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I. INTRODUCTION 



1. Introduction 

Since the experimental demonstration of tomographic quantum-state measurement on 
optical fields by Smithey, Beck, Raymer and Faridani [1993] based on the theoretical 
work of Vogel, K., and Risken [1989] a boom in studying quantum-state reconstruction 
problems has been observed. Numerous workers have considered various schemes and 
methods for extracting all knowable information on the state of a quantum system 
from measurable data, and early work has been recovered. In the history of quantum 
theory the concept of quantum state including quantum-state measurement has been 
a matter of intense discussion. 

The state of a quantum object is commonly described by a normalized Hilbert- 
space vector or, more generally, by a density operator g = J2i Pi\^i){^i\ which is a 
Hermitian and non-negative valued Hilbert-space operator of unit-trace. The Hilbert 
space of the object is usually spanned up by an orthonormalized set of basic vectors 
| A) representing the eigenvectors (eigenstates) of Hermitian operators A associated 
with a complete set of simultaneously measurable observables (physical quantities) of 
the object. 1 The eigenvalues A of these operators are the values of the observables 
that can be registered in a measurement. Here it must be distinguished between 
an individual (single) and an ensemble measurement (i.e., in principle, an infinitely 
large number of repeated measurements on identically prepared objects). Performing 
a single measurement on the object, a totally unpredictable value A is observed in 
general, and the state of the object has collapsed to the state \A) according to the 
von Neumann's projection definition of a measurement (von Neumann [1932]). If 
the same measurement is repeated immediately after the first measurement (on the 
same object), the result is now well predictable - the same value A as in the first 
measurement is observed. Obviously, owing to the first measurement the object has 
been prepared in the state \A). Repeating the measurement many times on identically 
prepared object, the relative rate at which the result A is observed approaches the 
diagonal density-matrix element as the number of measurements tend to 

infinity. Measuring for all values of A the statistics of A (and any function 

of A) are known. To completely describe the quantum state, i.e., to determine all the 
quantum-statistical properties of the object, knowledge of all density-matrix elements 
is needed. In particular, the off-diagonal elements essentially determine the 
statistics of such sets of observables B that are not compatible with A (\A,B\ ^ 
0) and cannot be measured simultaneously with A. Obviously, the statistics of B 
can also be obtained directly - similar to the statistics of A - from an ensemble 
measurement yielding the diagonal density- matrix elements (B\g\B) in the basis of 
the eigenvectors \B) of B. Now one can proceed to consider other sets of observables 
that are not compatible with A and B, and the rather old question arises of which 
and how many incompatible observables must be measured in order to obtain all the 
information contained in the density matrix of an arbitrary quantum state. Already 
Pauli [1933] addressed the question whether or not the wave function of a particle 
can be reconstructed from the position and momentum probability distributions, i.e., 
the diagonal density-matrix elements in the position and momentum basis. 

Roughly speaking, there have been two routes to collect measurable data for 

1 For notational convenience we write A, without further specifying the quantities belonging to 
the set. 
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reconstructing the quantum state of an object under consideration. In the first, 
which closely follows the line given above, a succession of (ensemble) measurements 
is made such that a set of noncommutative object observables is measured which 
carries the complete information about the quantum state. A typical example is the 
tomographic quantum-state reconstruction mentioned at the beginning. Here, the 
probability distributions p(x,ip) = (x,ip\g\x,ip) of the quadrature components x(ip) 
of a radiation- field mode for all phases ip within a n interval are measured, i.e., 
the expectation values of the quadrature-component projectors \x,ip)(x,ip\. In the 
second, the object is coupled to a reference system (whose quantum state is well 
known) such that measurement of observables of the composite system corresponds 
to "simultaneous" measurement of noncommutative observables of the object. In this 
case the number of observables (of the composite system) that must be measured in a 
succession of (ensemble) measurements can be reduced drastically, but at the expense 
of the image sharpness of the object. As a result of the additional noise introduced by 
the reference system only fuzzy measurements on the object can be performed, which 
just makes a "simultaneous" measurement of incompatible object observables feasible. 
A typical example is the Q function of a radiation-field mode, which is given by the 
diagonal density-matrix elements in the coherent-state basis, Q(a) = n^ 1 (a\g\a) . It 
can already be obtained from one ensemble measurement of the complex amplitude 
a, which corresponds to a fuzzy measurement of the "joint" probability distribution 
of two canonically conjugated observables x(ip) and x(ip + 7r/2) of the object. 

Since the sets of quantities measured via the one or the other route (or an appro- 
priate combination of them) carry the complete information on the quantum state 
of the object, they can be regarded, in a sense, as representations of the quantum 
state, which can be more or less close to (or far from) familiar quantum-state repre- 
sentations, such as the density matrix in the Fock basis or the Wigner function. In 
any case, the question arises of how to reconstruct from the measured data specific 
quantum-state representations (or specific quantum-statistical properties of the ob- 
ject for which a direct measurement scheme is not available). Again, there have been 
two typical concepts for solving the problem. In the first, equations that relate the 
measured quantities to the desired quantities are derived and tried to be solved either 
analytically or numerically in order to obtain the desired quantities in terms of the 
measured quantities. In practice, the measured data are often incomplete, i.e., not all 
quantities needed for a precise reconstruction are measured, 2 and moreover, the mea- 
sured data are inaccurate. Obviously, any experiment can only run for a finite time, 
which prevents one, on principle, from performing an infinite number of repeated 
measurements in order to obtain precise expectation values. These inadequacies give 
rise to systematic and statistical errors of the reconstructed quantities, which can 
be quantified in terms of confidence intervals. In the second, statistical methods are 
used from the very beginning in order to obtain the best a posteriori estimation of the 
desired quantities on the basis of the available (i.e., incomplete and/or inaccurate) 
data measured. However the price to pay may be high. Whereas in the first concept 
linear equations are typically to be handled and estimates of the desired quantities 
(including statistical errors) can often be directly sampled from the measured data, 
application of purely statistical methods, such as the principle of maximum entropy or 



2 To compensate for incomplete data, some a priori knowledge of the quantum state is needed. 
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Bayesian inference, require nonlinear equations to be considered, and reconstruction 
in real time is impossible in general. 

The aim of this review is to familiarize physicists with the recent progress in the 
field of quantum-state reconstruction and draw attention to important contributions 
in the large body of work. Although in much of what follows we consider the recon- 
struction of quantum states of travelling optical fields, the basic-theoretical concepts 
also apply to cavity fields and matter systems. Section 2. reviews the experimen- 
tal schemes that have been considered for quantum-state measurement on optical 
fields, with special emphasis on optical homodyning (§2.1.). It begins by formulat- 
ing the basic ideas of four-port homodyne detection of the quadrature components 
of single-mode fields (§§ 2.1.1. and 2.1.2.) and then proceeds to extend the scheme 
to multimode fields and multiport homodyne detectors (§ 2.1.3.), which can also be 
used for "simultaneous" measurement of noncommuting signal observables (§§ 2.1.4. 

- 2.1.6.). Section 2.1.7. returns to the four-port scheme in § 2.1.1. and explains its use 
for measuring displaced Fock states. Homodyne correlation measurements are shown 
to yield insight in phase-sensitive field properties even in the case of low detection 
efficiencies (§ 2.1.8.). After addressing heterodyne detection (§ 2.2.) and paramet- 
ric amplification (§ 2.3.), typical schemes for quantum-state measurement in high-Q 
cavity QED by test atoms are outlined (§ 2.4.). 

Section 3. reviews typical methods for reconstruction of quantum-state represen- 
tations and specific quantum-statistical properties of optical fields from sets of mea- 
surable quantities carrying the complete information on the quantum state. In § 3.1. 

- 3.5. we focus on the reconstruction of quantum-state representations from quantities 
measurable in homodyne detection of travelling optical fields. The next subsection 
(§ 3.6.) presents methods of reconstruction of quantum-state representations of high- 
er cavity fields from measurable properties of test atoms probing the cavity fields, and 
in § 3.7. alternative proposals are outlined. Section 3.8. then addresses the problem 
of direct reconstruction of specific quantities from the measured data - an important 
problem with respect to quantities for which direct measurement schemes have not 
been found so far. Whereas in §§ 3.1. - 3.8. it is assumed that all the data needed for 
a precise reconstruction of the desired quantities can be precisely measured, at least 
in tendency, in § 3.9. statistical methods for processing smeared and inaccurately 
measured incomplete data are outlined with the aim to obtain optimum estimations 
of the quantum states. 

Section 4. summarizes methods for measurement and reconstruction of quantum- 
state representations of matter system, with special emphasis on optical methods 
and experimental demonstrations. In § 4.1. the problem of reconstruction of the 
vibrational quantum state in an excited electronic state of a two-atom molecule is 
addressed. Typical schemes for determining the motional quantum state of trapped 
atoms are outlined in § 4.2., and in § 4.3. the problem of quantum-state measure- 
ment on Bose-Einstein condensates is considered. Sections 4.4. and 4.5., respectively, 
present schemes suitable for determining the motional quantum state of atom- and 
electron- wave packets. In the debate on quantum-state measurement, which is nearly 
as old as quantum mechanics, spin systems has been attracted much attention, be- 
cause of low-dimensional Hilbert space. Section 4.6. outlines the main ideas in this 
field. Finally, § 4.7. explains the basic ideas of Compton and X-ray scattering for 
measuring the single-particle quantum state of crystal lattices. 
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2. Phase-sensitive measurements of light 

The statistical properties of a classical radiation field are known when the amplitude 
statistics and the phase statistics are known. Whereas the amplitude statistics can 
be obtained from intensity measurements, determination of the phase statistics needs 
phase-sensitive measurement. Obviously the same is true in quantum optics. In 
order to obtain the complete information about the quantum state a radiation field 
is prepared in interferometric measurements that respond to amplitude and phase 
variables must be performed. 

2.1. Optical homodyning 

Homodyne detection has been a powerful method for measuring phase sensitive prop- 
erties of travelling optical fields which are suitable for quantum-state reconstruction, 
and a number of sophisticated detection schemes have been studied. In the four- 
port basic scheme, a signal field is combined, through a lossless beam splitter, with 
a highly stable reference field which has the same mid-frequency as the signal field. 
The reference field, also called local oscillator, is usually prepared in a coherent state 
of large photon number. The superimposed fields impinge on photodetectors, the 
numbers of the emitted (and electronically processed) photoelectrons being the ho- 
modyne detection output (for the basic ideas, see Yuen and Shapiro, J.H. [1978a,b, 
1980]; Shapiro, J.H., and Yuen [1979]; Shapiro, J.H., Yuen and Machado Mata [1979]). 
The observed interference fringes, which vary with the difference phase between the 
two fields, reflect the quantum statistics of the signal field and can be used - under 
certain circumstances - to obtain the quantum state of the signal field. The homo- 
dyne output can be fully given in terms of the joint-event probability distribution 
of the detectors in the output channels. In balanced homodyning, difference-event 
distributions are measured. In particular, the difference-event statistics measurable 
by a perfect four-port homodyne detector directly yields the quadrature-component 
statistics of the signal field, which has offered novel possibilities of quantum-state 
measurement. 

2.1.1. Basic scheme 

Let us start the analysis with the four-port scheme and first restrict attention to 
single-mode fields mode-matched to the local-oscillator (Fig. I). 3 The action of a 
lossless beam splitter (as an example of a linear four-port coupler) on quantum fields 
has widely been studied (see, e.g., Richter, G., Brunner and Paul, H., [1964]; Brunner, 
Paul, H., and Richter, G., [1965]; Paul, H., Brunner and Richter, G., [1966]; Yurke, 
McCall and Klauder [1986]; Prasad, Scully and Martienssen [1987]; Ou, Hong and 
Mandel [1987a]; Fearn and Loudon [1987]; Campos, Saleh and Teich [1989]; Leonhardt 
[1993]; Vogel, W., and Welsch [1994]; Luis and Sanchez-Soto [1995]), and it is well 
known that the input-output relations can be characterized by the SU(2) Lie algebra. 

In the Heisenberg picture, the photon destruction operators of the outgoing modes, 
a' k (k = 1,2), can be obtained from those of the incoming modes, dk, by a unitary 



3 Here and in the following spatial-temporal modes are considered. They are nonmonochromatic 
in general and defined with a particular spatial, temporal and spectral form (Appendix A.). 



8 



2. PHASE-SENSITIVE MEASUREMENTS OF LIGHT 





j dt 




ADC 








Figure 1: A signal pulse field E s interferes with a shorter local-oscillator pulse 
Ei at a 50/50 beam splitter. The local-oscillator phase <fi determines which 
quadrature amplitude of the signal is detected. The superposed fields are 
detected with high-efficiency photodiodes having response times much longer 
than the pulse durations. The photocurrents are integrated and sampled by 
analog-to-digital converters, to yield pulse photoelectron numbers n\ and n 2 . 
(After Raymer, Cooper, Carmichael, Beck and Smithey [1995].) 



transformation 



k'=l 



where Uu and f/ 2 i are the transmittance and reflectance of the beam splitter from one 
side, and U 22 and U 12 are those from the other side. The unitarity of the scattering 
matrix Ukk' — \Ukk'\^ kk ' implies that the following restrictions are placed on it: 

Wu\ 2 + \U 21 \ 2 = 1, \U 22 \ = |C/ n |, \U 12 \ = |C/ 21 |, (2) 

V?ll - V21 + V?22 - f!2 = ±7T. (3) 

The relations (2) and (3) ensure that the bosonic commutation relations for the pho- 
ton creation and destruction operators are preserved. Using the angular momentum 
representation, the unitary transformation (1) can be given by 

a' k = Va k V\ (4) 

where the operator V reads as 4 

V = e -^s e -ipL 2e -ijL 3e -i6N_ ^ 

Here, N = h\ + h 2 is the operator of the total photon number {hk = a\.ak), and the 
operators 

Li = 7; (a\a 2 + a\aij , L 2 = ^ (a\a 2 — a^i) , L 3 — \ (a\di — d 2 a 2 ) (6) 
obey the familiar commutation relations of angular-momentum operators. 



4 The parameters in eq. (5) are related to the parameters in eqs. (2) and (3) as a = \{<pu — ^22 
+ </?2i -v?i2±7r), 7 = i(^n-^22-^2i + <^i2T7r), cos 2 (0/2) = \ U n 1 2 , sin 2 (/3/2) = |[/ 2 i| 2 , 6= ±(<pn 
+ ^22)- 
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In the Schrodinger picture the photonic operators are left unchanged, but the 
state is transformed. Let g and g' be the input- and output-state density operator of 
the two modes, respectively. It is easily seen that if they are related to each other as 

g' = V^gV, (7) 

then the two pictures lead to identical expectation values of arbitrary field quantities. 

Let us now assume that the operation of the photodetectors in the two output 
channels of the beam splitter can be described by means of the standard photode- 
tection theory as outlined in Appendix C; i.e., the numbers of photoelectric events 
counted in a given time interval (integrated photocurrents) are proportional to the 
numbers of photons falling on the detectors. From eq. (C.8) the joint-event distribu- 
tion of measuring mi and m 2 events in the two output channels 1 and 2, respectively, 
is then given (in the Heisenberg picture) by 



where 




with d' k being given by eq. (1) [or eq. (4)]. Note that for rj — 1 (i.e., perfect detection) 
-fmi,m 2 1S nothing more than the joint-photon- number distribution of the two outgoing 
modes. In particular when the reference mode is prepared in a coherent state \oll), 
then in the expectation values of normally ordered operator functions, d' k can be 
replaced with (a = ai) 

a' k = U kl a + U k2 a L . (10) 



2.1.2. Quadrature-component statistics 

Let us now turn to the question of what information on the quantum statistics of the 
signal mode can be obtained from the homodyne output. From eq. (9) together with 
eqs. (10) and (3) we easily find that the output photon number n' k can be rewritten 

as 

n' k = \U k2 \ 2 \a L \ 2 + \U kl \ 2 tfa + (-l) k+1 \U k2 \\U kl \F(cp), (11) 

where 

F(<p) = 2 1 t 2 \a L \x{<p), (12) 

with 

x((p) = 2- 1/2 (a e'* + tie*) (13) 
being the quadrature-component operator, and 

(P = (Pl + <Pi2-<Pii- (14) 

We see that n' k consists of three terms. The first and the second terms represent the 
photon numbers of the local oscillator and the signal, respectively. More interestingly, 
the third term that arises from the interference of the two modes is proportional to 
a quadrature component of the signal mode, x(tp), the rapidly varying optical phase 
cut that usually occurs in a field strength (quadrature component) is replaced with 
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the phase parameter ip. In particular, when the signal and the local oscillator come 
from the same source, then the phase parameter if can be controlled easily, so that 
the dependence on <p of x(ip) can be monitored; e.g., by shifting the phase difference 
between the signal and the local oscillator in the input ports of the beam splitter in 
Fig. 1. 

If a balanced (i.e., 50%:50%) beam splitter is used, the relation |C4 2 | = |£4i|= l/\/2 
is valid, and hence the difference photon number is a signal- mode field strength F((p), 

h' 1 -fii l = F( ( p); (15) 

i.e., it is proportional to a signal-mode quadrature component. This result suggests 
that it is advantageous to use a balanced scheme and to measure the difference events 
or the corresponding difference of the photocurrents of the detectors in the two out- 
put channels in order to eliminate the intensities of the two input fields from the 
measured output. The method is also called balanced homodyning and can be used 
advantageously in order to suppress perturbing effects due to classical excess noise 
of the local oscillator (Mandel [1982]; Yuen and Chan [1983]; Abbas, Chan and Yee 
[1983]; Schumaker [1984]; Shapiro, J.H., [1985]; Yurke [1985]; Loudon [1986]; Collett, 
Loudon and Gardiner [1987]; Yurke and Stoler [1987]; Yurke, Grangier, Slusher and 
Potasek [1987]; Drummond [1989]; Blow, Loudon, Phoenix and Shepherd [1990]; Hut- 
tner, Baumberg, Ryan and Barnett [1992]). In particular, when the local oscillator is 
much stronger than the signal field, then even small intensity fluctuations of the lo- 
cal oscillator may significantly disturb the signal-mode quadrature components that 
are desired to be observed. For example, if the quadrature-component variance is 
intended to be derived from the variance of events measured by a single detector, the 
classical noise of the local oscillator and the quantum noise of the signal would con- 
tribute to the measured data in the same manner, so that the two effects are hardly 
distinguishable. Since in the two output channels identical classical-noise effects are 
observed, in the balanced scheme they eventually cancel in the measured signal due 
to the subtraction procedure. 

Strictly speaking, from the arguments given above it is only established that for 
chosen phase parameters the mean value of the measured difference events or pho- 
tocurrents is proportional to the expectation value of a quadrature component of the 
signal mode. From a more careful (quantum-mechanical) analysis it can be shown that 
in perfect balanced homodyning the quadrature-component statistics of the signal 
mode are indeed measured, provided that the local oscillator is sufficiently strong; i.e., 
\cil\ 2 is large compared with the mean number of signal photons (Carmichael [1987]; 
Braunstein [1990]; Vogel, W., and Grabow [1993]; Vogel, W., and Welsch [1994]; 
Raymer, Cooper, Carmichael, Beck and Smithey [1995]; Munroe, Boggavarapu, An- 
derson, Leonhardt and Raymer [1996]). The probability distribution of the difference 
events, 

Am = mi — m 2 , (16) 
can be derived from the joint-event distribution (8) as 

Examples of the dependence of the difference-event distribution P& m on the local- 
oscillator strength and the detection efficiency are shown in Fig. 2. In the limit of 



(17) 
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the local oscillator being sufficiently strong, so that |«l| 2 is large compared with the 
average number of photons in the signal mode, the difference-event distribution P& m 
can be given by (assuming that the two detectors have the same quantum efficiency 

V) 

1 ( Am \ 

with the phase tp being determined by the phase parameters of the apparatus and of 
the local oscillator according to eq. (14). In eq. (18), p(x,(p;rj) is a convolution of 
the quadrature-component distribution of the signal mode, p(x, <p) = (x, <p\g\x, ip) [cf. 
Appendix B.], with a Gaussian: 5 



p(x,cp;r]) = J dx'p(x',cp)p(x-x';r]), 



(19) 



p(x;rj) 



rj 



cxp 



r]x 



(20) 



/ 7r(l — 77) ^ \ 1 — r] / 
In particular, for perfect detection, rj — 1, the measured difference-event distribution 
o.e 
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Figure 2: Difference statistics P\ n for a single-photon state, quantum effi- 
ciency 7] = 0.75, and various photon numbers of the local oscillator: |«l| 2 = 
(a), 0.5 (b), 5 (c), 10 (d). The solid curves represent the smeared quadrature- 
component distribution p(x, ip; rj), eq. (19), which is a convolution of p(x, ip) 
with p(x;r)). (After Vogel, W., and Grabow [1993].) 



-Pato is identical with the difference-photon-number distribution. In this case p(x; rf) 
reduces to a 5 function; i.e., p(x,<p;r]= 1) =p(x,<p), and hence P& m is (apart from 
a scaling factor) exactly a quadrature-component distribution p(x, ip) of the signal 
mode. Note that when the local oscillator is sufficiently strong, then Am/(\/2\a L \) 
in eq. (18) is effectively continuous. In other words, single-photon resolution is not 

5 This result is obtained directly (in the limit of a strong local oscillator) from the basic equations 
(8), (17) of photocounting theory (Vogel, W., and Grabow [1993]), which include the effects of 
imperfect detection. 
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needed in order to measure the (continuous) quadrature-component statistics with 
high accuracy. In this case, highly efficient linear response photodiodes can be used, 
which do not discriminate between single photons, but nearly reach 100% quantum 
efficiency (Polzik, Carri and Kimble [1992]). 

The Gaussian p(x; if) [with dispersion a 2 = (1 — 77)/ (277)] obviously reflects the noise 
associated with nonperfect detection (77 < 1). In particular, for 50% quantum efficiency 
(i.e., 77 = 1/2) the function p(x; 77) is the (phase-independent) quadrature-component 
distribution of the vacuum (cr 2 = 1/2). Obviously, the measured quadrature-compo- 
nent distribution p(x, tp; 77) does not correspond to the quadrature component x(ip) of 
the signal mode, but it corresponds to the quadrature component x(ip; 77) of a super- 
position of the quadrature component of the signal, x(<p), and that of an additional 
(Gaussian) noise source, xjv(^), 



(21) 



so that x(ip; 77) is effectively a combined two-mode quadrature component. In particu- 
lar, eq. (21) reveals that the effect of nonperfect detection can be modelled, assuming 
a (virtual) beam splitter is placed in front of a perfect detector, since eq. (21) exactly 
corresponds to a beam-splitter transformation (§ 2.1.1.) (Yurke and Stoler [1987]; 
Leonhardt and Paul, H., [1993b, 1994a,b]). In this case the fields are only partly de- 
tected, together with some fraction of vacuum noise introduced through the "unused" 
input ports of the beam splitters. 6 

If the signal field is not perfectly mode-matched to the local-oscillator mode, the 
non-mode-matched part can give rise to additional noise in the measured quadrature- 
component distributions of the mode-matched signal (Raymer, Cooper, Carmichael, 
Beck and Smithey [1995]). It can be shown that when the local oscillator is strong 
enough to dominate the mode-matched signal field but not necessarily the non- 
matched signal, then the detection efficiency reduces to 



77 = 77 1 + 



N 



B 



\ a L\ 



(22) 



where N b is the mean number of photons in the nonmatched signal; i.e., 77 is replaced 
with 77 in eq. (18). As expected, 77 approaches 77 only in the limit \&l\ 2 Nb- 

In terms of the characteristic functions of the distributions involved, the content 
of eqs. (18) - (20) can be given by 



^a(v) = ^{V^y,f,v) = exp 



1)77 + 1^ 2 



277 



Here, 



(Am N 
i~, \V 



$(V^;s). (23) 



(24) 



is the characteristic function of the measured probability distribution of the scaled 
difference events Am/ (r}\aL\), and ^(z, ip; 77) and s), respectively, are the charac- 
teristic functions, i.e., the Fourier transforms, of the (owing to nonperfect detection) 



6 For a discussion of eqs. (19) and (20) in terms of moments of the quadrature components x(ip; rj) 
and x((p), see Banaszek and Wodkiewicz [1997a]. 
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smeared quadrature-component distribution p(x, ip; rj) and the s-parametrized phase- 
space function P(a; s) [see eqs. (B.26) 7 and (B.23)]. Equation (23) reveals that in the 
case of perfect detection the characteristic function of the (scaled) difference-event 
distribution is equal, along a line in the complex plane, to the characteristic function 
of the Wigner function, which can be regarded as a suitable representation of the 
quantum state of the signal mode. Hence, when the difference-event statistics are 
measured on a sufficiently dense grid of phases (p within a n interval, then all know- 
able information on the signal-mode quantum state can be obtained. The result is 
obviously a consequence of the fact that knowledge of the quadrature-component dis- 
tribution for all phases within a n interval is equivalent to knowledge of the quantum 
state; see Appendix B.5. (Vogel, K., and Risken [1989]). 

2.1.3. Multimode detection 

Relations of the type given in eq. (23) are also valid for (higher than four-port) 
multiport homodyning (Walker [1987]; Kiihn, Vogel, W., and Welsch [1995]; Ou and 
Kimble [1995]). Let us consider a linear, lossless 2(iV+ l)-port coupler (Fig. 3) and 
assume, e.g., that N spatially separated signal modes (channels 1, . . . , N) and a strong 
local-oscillator mode (channel N+l) of the same frequency are mixed to obtain (JV+1) 
output modes impinging on photodetectors. 8 The input and output photon operators 
a k and a' k , respectively, are then related to each other by a SU(iV-l-l) transformation, 
extending eq. (1) to 

JV+l 

a' k = U kk , a k , (25) 
k'=l 

(k — 1, . . . , N + 1). Note that any discrete finite-dimensional unitary matrix can be 
constructed in the laboratory using devices, such as beam splitters, phase shifters and 
mirrors (Reck, Zeilinger, Bernstein and Bertani [1994]; for further readings, see Jex, 
Stenholm and Zeilinger [1995]; Mattle, Michler, Weinfurter, Zeilinger and Zukowski 
[1995]; Stenholm [1995]; Torma, Stenholm and Jex [1995]; Torma, Jex and Stenholm 
[1996]). Simultaneous detection of the output modes yields an (JV + l)-fold joint 
event distribution P mi ,...,m N+1 - With regard to difference-event measurements, let us 
consider the scaled difference events 

Am, = 4VM^±4^. (26) 

m\ a L\ \u kN +i\ m\<XL\ 

Here, l — l,... ,k — l,k + l,. . . ,N + 1, the reference channel being denoted by k. It can 
be shown that when the photon number \ai\ 2 of the coherent-state local oscillator is 
much greater than unity and dominates all the signal modes at the detectors, 9 then 



7 The characteristic function \I> (z , ip; rf) of p(x,<p\rj) is denned according to eq. (B.26) with 
p(x, if; rj) in place of p(x, p). It is seen to be related to the characteristic function ^(z, ip) of p(x, ip) 
as p\ rf) =cxp[— i(?7 _1 - l)z 2 ] ^{z, p). 

8 All the modes are assumed to be matched to each other. 

9 Apart from this condition, the analysis also allows the local oscillator to be in other than coherent 
states (Walker [1987]; Leonhardt [1995]; Kim, M.S., and Sanders [1996]; for measuring the degree 
of squeezing in a signal field by using squeezed local oscillators, see Kim, C, and Kumar [1994]). 
This can be done by writing the local-oscillator state ql in the form ql = D{ai) q' l & {a^) , where 
D{oil) is the coherent displacement operator and g' L is a state that remains finite as the complex 
amplitude |aj,| becomes large. 
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Figure 3: Scheme of homodyne 2(iV + l)-port detection. N signal modes 
a k (k — l,...,N) of frequency uj and a strong local oscillator cll of the same 
frequency are combined by a lossless 2(iV + l)-port device to give N + 1 output 
modes (U kk >: unitary transformation matrix of the device). Simultaneous 
detection of the N+1 output modes yields the (N + l)-fold joint count 



distribution P, 



mi, 



,m N - 



(After Vogel and Welsch [1995].) 



the characteristic function of the joint probability distribution of the scaled difference 
events, £Ia{{vi}), can be related to the characteristic function of the s-parametrized 
multimode phase-space function of the signal field, $({/9 n }; s), n — 1, . . . , N, as 10 



^a({vi}) = exp 



(s-l)77+l 



N 



EIA.I 



<&({/?„}; s), (27) 



2^ n=1 

where equal detection efficiencies have been assumed (i] k = i] k r =77), and 

N+1 N+1 

Pn = ie ilfiL E U kN+1 U* kn y k , with ^ \U kN+l \ 2 y k = (28) 

k=l k=l 

(Kiihn, Vogel, W., and Welsch [1995]). Note that when N — l (four-port apparatus) 
and \U k 2 \ = \U k i \ = l/\/2, then eq. (27) [together with eq. (28)] reduces to eq. (23). 
It should be pointed out that joint measurements on combinations of multiports of 
the above described type can also be used in order to detect (groups of) correlated 
modes of different frequencies. In particular, the combination of two balanced four- 
port schemes can be used to measure the joint quadrature-component distribution 
p(x 1 ,x 2 , (fii, V2) °f a (correlated) two- mode signal field for all values of the two phases 
(pi and (pi within 7r intervals in order to determine the two-mode quantum state 
(Raymer, Smithey, Beck, Anderson and McAlister [1993]). 



10 Note that $({/?,}; s) = exp \ s |/3„| 2 ] ( D{{(3 n }) ), where D{{f3 n }) is the iV-mode coherent 
displacement operator, and the Fourier transform of $({/?„}; s) is the s-parametrized iV-mode phase- 
space function. 
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From eqs. (27) and (28) it can be seen that for perfect detection the character- 
istic function of the iV-fold joint scaled difference-event distribution is nothing other 
than the characteristic function of the Wigner function of the iV-mode input radi- 
ation field at certain values of the N complex arguments /3 n , and it corresponds to 
the characteristic function of a joint quadrature-component distribution. To obtain 
all knowable information on the iV-mode quantum state, i.e., the joint quadrature- 
component distributions for all relevant phases or the complete Wigner function, the 
arguments j3 n must be allowed to attain arbitrary complex values. This may be 
achieved by means of an appropriate succession of (ensemble) measurements or one 
(ensemble) measurement, including in the apparatus appropriately chosen reference 
modes whose quantum states are known (or a combination of the two methods) (Vo- 
gel, W., and Welsch [1995]). In the first case, phase shifters in the apparatus may 
be used to appropriately vary (similar to the four-port scheme) the scattering matrix 
Ukk' from measurement to measurement. In the second case, the scattering matrix is 
left unchanged, but additional reference inputs are used, so that a joint quadrature- 
component distribution of signal and reference modes is effectively measured (see 



Instead of measuring the iV-fold joint difference-event statistics in a 2(N + l)-port 
scheme, the information on the iV-mode signal field can also be obtained from the 
difference-event statistics measured in a standard four-port scheme with controlled 
signal-mode superposition (Opatrny, Welsch and Vogel, W., [1996, 1997a,b]; Raymer, 
McAlister and Leonhardt [1996]; McAlister and Raymer [1997a,b]). In the scheme, 
which applies to both spatially separable and nonseparable modes, the signal input 
is formed as a superposition of all the signal modes. Controlling the expansion coeffi- 
cients in the superposition, from the measured difference-event distributions the same 
information on the quantum statistics of the iV-mode signal-field can be obtained 
as from the iV-fold joint difference-event distribution measured in a 2(N + l)-port 
scheme. 

If the signal modes are separated spatially, they can be superimposed using an 
optical multiport interferometer such that each of the signal modes is fed into a 
separate input channel of the interferometer and the superimposed mode in one of 
the output channels is used as the (signal) input of a four-port homodyne detector. 
If the signal field is a pulse-like mode, the signal pulse can be combined, through 
the beam splitter in the four-port scheme, with a sequence of short local-oscillator 
pulses in order to analyse the signal pulse in terms of the modes associated with the 
local-oscillator pulses. In this case the superposition of the component modes to be 
measured is achieved by a proper formation of the local oscillator pulses. 

To illustrate the method, let us consider the simplest case when a (correlated) two- 
mode field is perfectly detected (Fig. 4). Then the measured quantitity is a weighted 
sum of the quadrature components Xi(<^i) and x 2 (y?2) of the two modes, 



where the phases ipi and y?2 and the relative weight [superposition parameter a, 
a G (0,7r/2)] are assumed to be controlled in the measurement. For spatially sep- 
arated modes [Fig. 4(a)] the phases can be controlled by phase shifters placed in 
front of a beam splitter used for mixing the modes, and the parameter a can be 



§ 2.1.4.). 




(29) 
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Figure 4: Two possible schemes for reconstructing two-mode density matrices 
via measurement of combined quadrature components, (a) Two modes (Si 
and a 2 ) are mixed by a beam splitter and one of the interfering output modes 
(6) is used as signal mode in balanced homodyning in order to measure the 
sum quadratures of the two modes (LO is the strong local oscillator), (b) A 
signal pulse and a sequence of two short (strong) local-oscillator pulses with 
envelopes Fi(t) and F 2 (t) are superimposed in balanced homodyne detection. 
(After Opatrny, Welsch and Vogel, W., [1997b].) 



controlled by the transmittance (reflectance) of the beam splitter. If the quadra- 
ture statistics of a signal pulse are measured in a set of short local-oscillator pulses 
[Fig. 4(b)], then the parameter a can be controlled by changing the mutual intensi- 
ties of the two pulses. Again, the phases can be controlled interferometrically. 11 Let 
Pj(x±,X2, (fii, (£2) and ps(x, a, ipi, ^2) be respectively, the joint quadrature-component 
distribution and the (measured) sum quadrature-component distribution of the two 
modes, and \& j(z±, Z2, <pi, V2) and &s(z, a , <Pi, V2) the corresponding characteristic 
functions. It can be shown that 

^ j(z cos a, z sin a, ip u <p 2 ) = ^ s (z,a,(p 1 ,(p 2 ). (30) 

In other words, the joint- quadrature component distribution for all phase values 
(within 7r intervals) can be obtained from the sum quadrature component distri- 
bution for all phase values and all values of the weighting factor, which shows the 
possibility to obtain the two-mode quantum state from the measured statistics of the 
summed quadratures. The effect of nonperfect detection can be taken into account 
according to eq. (23). 

2.1.4. Q function 

As mentioned in § 2.1.3., the 2A^-fold manifold of data necessary for a reconstruction of 
the quantum state of an iV-mode signal field can be obtained by including additional 
reference modes in the detection scheme in Fig. 3 whose quantum states are known, 
such as the vacuum inputs in "unused" input channels. Let us suppose that the 
2(N + l)-port apparatus is extended to a 2(2A^ + l)-port device, N input channels 

11 If signal and local oscillator come from different sources, only the phase difference 8<p= pi — f2 
can be controlled. 
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being "unused" . Since the vacuum inputs are not correlated to each other and to the 
signal input, the characteristic function of the quantum state of the (2iV)-mode input 
field can be factorized, 

2N 

*({&};*) = *({#}; «) *({&}; *), with *({&};*) = U J-vW 2 (3i) 

j=N+l 

(i — l,...,N, signal channels; j = N + 1, . . . , 27V, vacuum channels). Equation (28) 
(with 2N in place of N) can then be inverted easily in order to obtain 2N real 
arguments y\ of the characteristic function {{yi}) of the measured distribution for 
any (freely chosen) N complex arguments of the characteristic function <3>({/3j}; s) 
of the iV-mode signal-quantum state. In other words, in the extended detection 
scheme there is a one-to-one correspondence between the measured (2iV)-fold joint 
difference-event distribution and the quantum state of the iV-mode signal field (Vogel, 
W., and Welsch [1995]). 

Obviously, one local oscillator per group of (equal-frequency) signal modes and 
one additional reference input per signal mode are sufficient for measuring a joint 
difference-event distribution that can be regarded as a measure of the signal-mode 
quantum state. Needless to say, both the number of local oscillators and/or the 
number of additional reference inputs may be increased, and the quantum state of the 
signal may of course be expressed in terms of an appropriately chosen joint difference- 
event distribution recorded in such an extended measurement scheme. In the case of 
a single-mode signal field a six-port scheme is the minimum in order to determine the 
quantum state of the signal (Zucchetti, Vogel, W., and Welsch [1996]; Paris, Chizhov 
and Steuernagel [1997]). The situation is quite similar to classical optics, in which 
the six-port scheme is the minimum in order to determine the complex amplitude 
(Walker [1987]). Let us consider a 3 x 3 coupler and assume that the reference mode 
is in the vacuum state (Fig. 5). In this case, it can be shown that application of 
eqs. (27), (28) and (31) yields 

$(/3; s) = exp(| S |/?| 2 ) = exp{±|/?| 2 [s + (2-r^)^ 1 ] } Q A ( yi , y 2 ), (32) 

where (5 is relatated to y± and y 2 as 

P = li[(<f>*-l) yi + ((j>-l)y 2 ]e^, (33) 

with = exp(-i27r/3) [$(/3) = $(/3;0)]. Equation (33) reveals that $(/3; s = 1 - 2r]- r ) 
= Qa(ui, 2/2); i-e., measurement of the joint difference-event distribution ?Ami,Am 2 i n 
balanced six-port homodyning is equivalent to measurement of a phase-space function 
P(a; s) of the signal mode (in non-orthogonal coordinates in the phase space): 



\/3 
2r\ \OiL\ 



with 

a = -*?±Jj-^Jp (35) 

ands=l — 2i]~ 1 [y/<f) = exp(— iir/3)]. The function P(a; s = 1 — 2r]~ l ) can be regarded 
as a smoothed Q function of the signal mode, which approaches the Q function as the 
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unused port 
2^ 



signal field 
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Figure 5: Scheme of balanced six-port homodyning for the detection of the Q 
function. Three input fields (signal, local oscillator, vacuum) are combined by 
three symmetric beam splitters BSj (i — 1,2,3) and a (— 7r/2) phase shifter, 
where two (1/2 : 1/2) beam splitters (BSi and BS 3 ) and one (2/3 : 1/3) 
beam splitter (BS 2 ) are used. The joint difference statistics is recorded by the 
detectors PDj in the output channels [p s dc(-Di3, -D23) corresponds to ?A mi ,Am 2 
in eq. (34)]. (After Zucchetti, Vogel, W., and Welsch [1996].) 



quantum efficiency 77 goes to unity. In particular for perfect detection (77 = 1) the Q 
function of the signal mode is measured (in non-orthogonal phase-space coordinates). 

Identifying in eq. (B.25) s' withl — 2i]~ 1 , this equation can be regarded as a pre- 
scription for obtaining other phase-space functions P(a; s) from the measured one. 
When s < 1 — 2i]~ 1 is valid, then the j3 integration in eq. (B.25) can be performed 
separately to obtain P(a; s) as a convolution of P(a; s = 1 — 2r]~ 1 ) with a Gaussian, 
which reveals that all the phase-space functions which are typically broader than the 
measured one can be obtained simply by convolving the measured distribution with 
a Gaussian. In the opposite case, when s> 1 — 2r]~ 1 , then the integration over 7 must 
be done first in eq. (B.25). However, experimental inaccuracies may be exploding 
via the inverse Gaussian and prevent a stable reconstruction of P(a; s) with reason- 
able precision. Since the maximum value of s for an always stable reconstruction, 
i.e., s — 1 — 2r)~ 1 , tends to minus unity as 77 goes to unity, the upper boundary of s 
corresponds to the Q function. 

As mentioned above, higher than six-port homodyne detectors can also be used in 
order to measure the signal-mode quantum state in the phase space. Among them, 
eight-port homodyne detection has widely been studied (Walker and Carroll [1984]; 
Walker [1987]; Lai and Haus [1989]; Noh, Fougeres and Mandel [1991, 1992a,b]; Hradil 
[1992]; Freyberger and Schleich [1993]; Freyberger, Vogel, K., and Schleich [1993a,b]; 
Leonhardt and Paul, H., [1993a]; Luis and Pefina [1996a]). Let us consider two modes 
that are superimposed, in accordance with eq. (1), by a 50%:50% beam splitter, and 
assume that one of the incoming modes, say, the second is in the vacuum. It can 
then be shown that the joint probability distribution of the two outgoing quadrature 
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components x[((p) and x' 2 ((p + 7i/2) is the (scaled) Q function of the first incoming 
mode, 12 



p (x[, x' 2 , (p, v9 + 7r/2) ~ Q ai = 2 1/2 (x' 1 +ix 2 ) 



(36) 



(Lai and Haus [1989]; Leonhardt and Paul, H., [1993a]; for details, also see Leonhardt 
and Paul, H., [1995]; Leonhardt [1997c]). Hence, using the two outgoing modes 
of the beam splitter as incoming signal modes in two separate balanced four-port 
homodyne detectors and measuring (for a n/2 phase difference) the joint difference- 
event statistics of the two homodyne detectors, the Q function of the original signal 
mode is obtained, provided that perfect detection is accomplished. 13 Altogether the 
setup forms an balanced eight-port homodyne detector, with a signal input, a vacuum 
input and two local-oscillator inputs. Alternatively, the beam splitter and the two 
four-port homodyne detectors can be combined into an eight-port apparatus with two 
vacuum inputs and one local-oscillator input (Fig. 6). A straightforward calculation 



D 5 



□ A/4 



Figure 6: Balanced homodyne eight-port detection scheme for measuring the 
Q function, using 50%:50% beam splitters and a A/4 phase shifter. The 
signal is fed into port 1 and the strong local oscillator is fed into port 2, and 
vacuum is fed into the ports 1' and 2'. The joint difference-event probabbility 
distribution i-Ami,Am 2 > ec L- (37), is measured in the channels 3 and 4 (Ami) 
and 5 and 6 (Am 2 ). (After Vogel, W., and Welsch [1994].) 



12 Equation (36) can be proved correct, applying the beam splitter transformation (1) and ex- 
pressing the characteristic function of the two-mode (outgoing) quadrature-component distribution 
p (x[,x 2 , if, ip+n/2) in terms of the characteristic function of the single-mode (incoming) Q function 
Q(ai). 

13 Note that x(ip) and x((p + n/2) play the same role as position and momentum of a harmonic 
oscillator in quantum mechanics, because of [x(ip),x((p + ir/2)] = i. Therefore also the notations 
q = x(ip) and p = x(ip+Tr/2) are frequently used. Equation (36) is an example of "simultaneous" mea- 
surement of a pair of conjugate quantities. Actually, the quantities q[ = x'^ip) and p' 2 = x 2 (ip+ir/2) 
are measured, which can be regarded as the conjugate quantities q\ = Xi((p) and pi =xi(<p+n/2), 
respectively, "smoothed" by the introduction of additional (vacuum) noise necessary for a simul- 
taneous (but approximate) measurement of q\ and p\. Accordingly, an additional uncertainty is 
introduced in the measurement, and it can be shown that the uncertainty product for q[ and p' 2 
is twice that of the measurements of qi and p\ made individually, Aq[Ap 2 > 1 (Arthurs and Kelly 
[1965]; for a review, see Stenholm [1992]). 
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shows (similar to the six-port homodyne detector) that measurement of the joint 
difference-event distribution P/± mi Am 2 again yields the phase-space function P(a; s — 
1-2^), 

d 1 r>( Am 1 -iAm 2 A , 

"Ami,Am2 = n P \ a = J ] S = I - 2r] (37) 

(Freyberger and Schleich [1993]; Freyberger, Vogel, K., and Schleich [1993a,b]; Leon- 
hardt and Paul, H., [1993b]; Vogel, W., and Welsch [1994]; D'Ariano, Macchiavello 
and Paris [1995]; Kochahski and Wodkiewicz [1997]). Compared with the six-port 
scheme, PAmi,Am 2 is the (scaled) function P(a; s = 1 — 2rj~ 1 ) in an orthogonal ba- 
sis [cf. eqs. (34) and (37)]. The balanced eight-port homodyne detector was first 
used by Walker and Carroll in order to demonstrate the feasibility of measuring the 
components of the complex amplitude of optical signals, extending earlier microwave 
techniques (Walker and Carroll [1984]). 

Finally, it was proposed to measure the Q function by projection synthesis, mixing 
the signal mode with a reference mode that is prepared in a quantum state such 
that, for appropriately chosen parameters, the joint-photon-number probabilities in 
the two output channels of the beam splitter realize the coherent-state projector 
n(a) = ii~ l \a)(a\ for truncated signal states (Baseia, Moussa and Bagnato [1997]). 
The method was first introduced to synthesize (for truncated states) the phase-state 
projector fl(0) = (Barnett and Pegg [1996]; cf. § 3.8.3.). 



2.1.5. Probability operator measures 

As mentioned, the reference mode with which a signal mode is mixed has not neces- 
sarily to be in the vacuum state in order to obtain, in principle, all knowable infor- 
mation on the quantum state of the signal mode. If the reference mode is allowed 
to be prepared in a quantum state qr, then a joint measurement of the (7r/2-shifted) 
quadrature-components of the interfering fields can be regarded as a realization of 
a complex amplitude measurement (Walker [1987]). Each g R implies a probability 
operator measure (POM) over the complex amplitude which can be characterized by 
a positive valued Hermitian operator 

f[(a) = n^Dia) g R f)\a), (38) 

with 

Jd 2 afL(a) = I. (39) 

The joint probability density prob(a) of obtaining a result a = Re a + ilma from a 
measurement described by this POM is 14 

prob(a) = (f[(a)) . (40) 

14 The concept of POM was introduced in quantum theory in order to generalize the familiar 
quantum probability theory based on orthogonal projectors (Davies, E.B., and Lewis [1970]; Holevo 
[1973]; Davies, E.B., [1976]; Helstrom [1976]). In particular, it can be used for describing "simul- 
taneous" measurement of noncommuting observables of an object, the corresponding POM being 
related to a quantum measurement of commuting observables in some extended Hilbert space. Note 
that for nonperfect detection (77 < 1) the photocounting formula (C.8) is also a POM. 
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Note that from the properties of 11(a) it follows that prob(a) > and / d 2 a prob(a) 
= 1. The operational probability density distribution prob(o;) in eq. (40), which is 
also called propensity (Popper [1982]), can be given by a convolution of the Wigner 
function W(a) of the signal prepared in a state g with the Wigner function Wr(o) 
of the reference mode prepared in a state qr, 

prob(a) = J d 2 (3W R ((3-a)W((3) (41) 

(Husimi [1940]; Arthurs and Kelly, Jr., [1965]; Kano [1965]; Wootters and Zurek 
[1979]; O'Connell and Rajagopal [1982]; Rajagopal, A.K., [1983]; Wodkiewicz [1984, 
1986, 1988]; Takahashi and Saito [1985]; Walker [1987]; Lai and Haus [1989]; Hradil 
[1992]; Lalovic, Davidovic and Bijedic [1992]; Stenholm [1992]; Davidovic, Lalovic and 
Bijedic [1993]; Leonhardt [1993]; Leonhardt and Paul, H., [1993b, 1995]; Chaturvedi, 
Agarwal and Srinivasan [1994]; Raymer [1994]; Buzek, Keitel and Knight [1995a,b]; 
Paris, Chizhov and Steuernagel [1997]; Wunsche and Buzek [1997]). Obviously, the 
reference mode acts as a filter and smoothes the Wigner function of the signal. A filter 
system of this type is also called quantum ruler (Aharonov, Albert and Au [1981]). 
It is needed in order to resolve the noncommutative quadrature components of the 
signal. 

The particular realization of the filter strongly influences the outcome of the mea- 
surement. The class of phase-space functions that can be obtained includes the s- 
parametrized functions with s< — 1. In particular, detection of the Q-function implies 
that Qr = |0)(0|, so that 11(a) = 7r _1 |o;)(a|, and hence prob(o;) = 7r^ 1 (o;|^|a) = Q(a). 
This is the case when vacuum is fed into the reference channels in the homodyne 
detection schemes in Figs. 5 and 6, and the signal is mixed with it. If the signal is 
mixed with a squeezed vacuum, then the POMs of the form given in eq. (38) with 
^ R = 1 S(^)|0)(0| 1 S t (O can be realized (Walker [1987]; Lai and Haus [1989]; Leonhardt 
[1993]; Kim, M.S., and Sanders [1996]; Kochahski and Wodkiewicz [1997]). 15 

Probability operator measures over the complex amplitude and the photon number 
can also be defined such that the joint probability density for an outcome a = Reo;-|- 
ilm a and n is given by 

prob(a, n) = (fl(a, n)^ , (42) 

where Il(a,n) must be a positive valued Hermitian operator that can be used to 
resolve the idendity, 

J2 I d 2 all(a,n) = /. (43) 

A probability operator measure of this form can be realized by ten-port homodyne 
detection (Luis and Pefina [1996b]). 16 In the scheme, a signal mode is fed into one 
of the input channels of an unbalanced beam splitter whose second input channel is 
"unused". The mode in one of the output channels of the beam splitter (say, the 
transmitted signal) is used as the input mode in an eight-port homodyne detector, 

15 It should be noted that these POMs can also be realized in heterodyne detection (§ 2.2.) (Yuen 
and Shapiro, J.H., [1980]; Yuen [1982]). Further, they can be realized in an unbalanced homodyne 
scheme with vacuum input but not equal-part signal-beam splitting (Leonhardt [1993, 1997c]). 

16 The ten-port scheme can be regarded as a reduced twelve-port scheme, the latter being suited 
for measuring simultaneously the Stokes parameters of a two-mode field. 
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and a photodetector is placed in the other output channel in order to measure the 
photon- number statistics of the mode in this channel, 17 together with the eight-port 
complex-amplitude measurement. The realized POM is related closely to photon- 
added coherent states, 

U(a,n) = g n (a) a) n \a)(a\ a n , (44) 

where g n (oi) = 7r _1 t -2 exp[|a| 2 (l — t~ 2 )] (1 — t 2 ) n /n\ , with t being the (absolute value 
of the) transmittance of the beam splitter (0 < t < 1). Note that when t — > 1 then 
tl(a,n) — > 7T — 1 |ct) {ct| S n0 ; i.e., the familiar coherent-state POM is realized. 

2.1.6. Positive P function 

It is well known that for s > — 1 the s-parametrized phase-space functions P(a; s) 
(see Appendix B.4.) do not necessarily exist as positive functions, and for s > they 
are not necessarily well-behaved. The latter also applies to the P function, which is 
widely used to calculate averages of normally ordered (i.e., measurable) quantities. 
In order to avoid using singular functions, generalized P representations may be used 
(Drummond and Gardiner [1980]; see also Gardiner [1983, 1991]), such as the positive 
P function, 18 

P(a, a') = ^ exp(-i|o! - a'*\ 2 ) (±(a + a'*)\ g\±(a + a'*)). (45) 

The possibility of measuring the single-mode Q function in perfect (six-port or 
eight-port) homodyne detection offers also the possibility of measuring the quantum 
state of the mode in terms of the positive P distribution using more involved multiport 
homodyne detection schemes (Agarwal and Chaturvedi [1994]). Let us again consider 
a signal mode and a reference mode that are superimposed by a 50%:50% beam 
splitter and assume that the reference mode is in the vacuum. It can then be shown 
that the joint Q function of the two outgoing modes, Q (a±, 0:2), is related to the Q 
function of the signal mode, Q(cx), as 19 



Q(a 1 ,a 2 ) = exp(-i|a:i-ai2| 2 ) Q^a = 



(46) 



which reveals that Q(ai,a 2 ) is nothing but the (scaled) positive P function of the 
signal mode, 

Q(a 1 ,a 2 )=4P(a = V2a 1 , (3 = V2a* 2 ). (47) 

Hence, if each of the two output modes of the beam splitter is used as an input mode 
of a multiport apparatus (such as the six- or eight-port homodyne detector outlined 
in § 2.1.4.) that measures the Q function, then measurement of the joint Q function of 
the two modes yields the positive P function of the signal-mode under study. Needless 
to say that for imperfect detection a smeared positive P function is measured. 



17 For a direct measurement of the photon-number statistics, see § 2.1.7.. 

18 Note that g= J d 2 a J d 2 a' \a)(a'*\ ((a'*|a)) _1 P(a,a') in this representation. 

19 Equation (46) can be proved correct, applying the beam splitter transformation (1) and ex- 
pressing the characteristic function of the two- mode (outgoing) Q function Q(ai, 0:2) in terms of the 
characteristic function of the signal-mode (incoming) Q function Q(a). Note that Q(ai,a2) is the 
two-mode s-parametrized phase-space function P(ai,a 2 ; s = —1). 
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The positive P function is an example of a measurable phase-space function 
that is defined as a function of two complex amplitudes a and (5 (per mode). The 
concept can also be extended - similar to § 2.1.5. - to other than vacuum reference 
inputs (in the beam splitter and/or the multiport homodyne detectors) in order to 
obtain generalized phase-space functions defined in a four-dimensional phase space. 
A simple example is the smeared positive P function mentioned. 



2.1.7. Displaced-photon-number statistics 

Let us return to the four-port homodyne detection scheme (Fig. 1) and answer the 
question of which quantity is measured when a signal mode is mixed with a local- 
oscillator mode by an unbalanced beam splitter and only a single-channel homodyne 
output is measured. We rewrite eq. (10) as 

a' k = U kl (a -a) = U kl D(a)af)\a), (48) 

where 

a = -—a L , (49) 
Ukl 

and apply the photocounting formula (C.5). The probability of detecting m events 
in the kth output channel can then be given by 



P m = (:—Mn(a)] m e-^ :\ , (50) 



ml 



where n(a) = D{a)hD\a) is the displaced photon-number operator of the signal 
mode, and rj = \U k i\ 2 r)£) (with rjo being the quantum efficieny of the detector used). 
Equation (50) reveals that the observed probability distribution P m is nothing but 
the displaced photon-number distribution of the signal field measured with quan- 
tum efficiency rj. In particular, when a signal and a strong local oscillator, \at L \ 
— > oo, are mixed by a beam splitter with high transmittance, \Uu\ = | C/22 1 —> 1) and 
low reflectance, \U2i\ — \U±2\ — > 0, such that the product |£/i2||o!l| is finite, then for 
high quantum efficiency (77— > 1) the displaced photon-number probability distribu- 
tion of the signal is measured (Wallentowitz and Vogel, W., [1996a]; Banaszek and 
Wodkiewicz [1996]; Paris [1996a]), 

P m -> Pm{oi) = (m, a\g\m, a), (51) 

where \m,a) = D(a)\m) are the displaced photon-number states. It should be noted 
that for chosen m the quantity p m (a) as a function of a can be regarded (apart from 
the factor 7r _1 ) as a propensity for the signal-mode complex amplitude, which can 
also be measured in multiport homodyning (see § 2.1.5.). For chosen a it is an ordi- 
nary probability for the displaced signal-mode photon number, which can already be 
obtained from the four-port detector outlined here. In order to obtain in this scheme 
p m (a) as a function of a, a succession of (ensemble) measurements must be per- 
formed. Hence, measurement of the displaced-signal-mode photon-number statistics 



20 For a method suggested to measure the positive P function of a quantum-mechanical particle, 
see Braunstein, Caves and Milburn [1991]. 
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as a function of the complex parameter a is equivalent to measurement of the signal- 
mode quantum state, and it is expected that it yields more data than the minimum 
necessary for reconstructing it (§ 3.). 

In contrast to balanced homodyning, measurement of the displaced photon-num- 
ber statistics in unbalanced homodyning requires highly efficient photodetectors that 
can discriminate between n and n + 1 photons in order to resolve the discrete nature 
of the photon number. Presently, such detectors are not available. Photomultipliers 
and streak cameras can discriminate between single photons, provided that the field 
does not contain more than about 10 photons, but the quantum efficiency of about 
10% - 20% is extremely low. Currently available avalanche photodiodes operating in 
the Geiger regime may reach about 80% quantum efficiency (Kwiat, Steinberg, Chiao, 
Eberhard and Petroff [1993]), but they do not discriminate between single photons. 
They can only indicate the presence of photons, because of saturation. 

The problem may be overcome using multichannel coincidence-event measurement 
techniques also called photon chopping. In particular it was proposed to use highly 
efficient avalanche photodiodes and a beam splitter array to divide the number of 
readout photons among the photodiodes, so that none is likely to receive more than 
one photon (Ho, Lane, La Porta, Slusher and Yurke [1990]; Song, Caves and Yurke 
[1990]; Paul, H., Torma, Kiss and Jex [1996a]). 21 Alternatively, it was proposed to 
directly defocuse the field to be measured onto an array of photodiodes (Wallentowitz 
and Vogel, W., [1996a]). 

Let us assume that the mode to be detected enters one of the input ports of 
a balanced linear 2iV-port apparatus (the other input ports being "unused"), and 
multiple coincidences are measured at the output, placing avalanche photodiodes in 
the ./V output channels. If the signal mode contains less than N+l photons, then there 
is a one-to-one correspondence between the measured coincidence-event distribution 
p n (N) and the photon-number distribution p n of the signal. To be more specific, it 
can be shown that 22 

N 

Prn(N) = P m \n(N) Pn , (52) 

n=m 

where 

P H n(A0 = ^P|: (-iy(7)(™-*) n if rn<n<N, (53) 

and P m \ n (N) =0 if n < m (Paul, H., Torma, Kiss and Jex [1996a]). Here, P m \ n (N) 
is the probability of registration of m clicks under the condition that n photons are 
present. The probabilities P m \ n {N) form an N x N upper triangular matrix A mn = 
Pm\n(N) which can be inverted in order to calculate p n from p n (N), 

N 

Pn= J2 (^ _1 )nmPm- (54) 
m=n 



21 For detection of squeezing via coincidence-event measurement, see Janszky, Adam, P., and 
Yushin [1992]; for Fock state detection and preparation, see Paul, H., Torma, Kiss and Jex [1996b]. 
22 Here it is assumed that the balanced 2iV-port realizes a unitary transformation U^v = U 2 (8)11^/2 • 
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2.1.8. Homodyne correlation measurements 

As already mentioned, the accuracy with which a (quantum) field can be measured 
by a homodyne detector is limited by the overall quantum efficiency of the device [see, 
e.g., Eqs (19) and (20)]. To overcome this limitation, one may perform homodyne 
correlation measurements (Ou, Hong and Mandel [1987b]). In contrast to ordinary 
balanced homodyning, the (time-delayed) intensity correlation between the two out- 
going fields is measured. In particular, the information on squeezing is obtained from 
the time dependence of the measured correlation function. Since the measured coin- 
cidence events are proportional to the product of the detection efficiencies of the two 
detectors, small quantum efficiencies may reduce the measured signal (which could 
be compensated by longer measurement times) but do not smooth out the desired 
information. However, a drawback of the method is that the classical noise of the 
local oscillator is not balanced out, i.e., even small relative fluctuations of the (strong) 
local oscillator may prevent the quantum noise effects of a weak signal from being 
measured. 

It was therefore proposed to use a weak local oscillator whose intensity is com- 
parable to that of the signal (Vogel, W., [1991, 1995]). In this case, the classical 
noise of the (highly stablized) local oscillator may be reduced below the level of the 
quantum fluctuations of the signal. Moreover, simultaneous measurement of different 
kinds of correlation functions of the signal field is possible. To illustrate this, let 
us consider the difference between the measured second-order intensity correlation 
function G^(t, t + r) for short and long delay times r, 



and restrict attention to stationary fields, so that the time argument t can be omitted. 
Decomposing AG^ with respect to powers of the local-oscillator amplitude E L) one 
may observe the following effects. The zeroth-order term yields the normally ordered 
intensity (/) fluctuation of the signal, 



The second-order term is related to the normally ordered electric-field variance of the 
signal, 



where E(ip) corresponds to the quadrature- component operator x((p), <p being the 
phase difference between signal and local oscillator. Eventually, the first-order term 
represents the correlation between the two signal-field observables, 



AG (2) (t) = G {2) (t, t) - lim G {2) (t, t + r), 



(55) 




(56) 




(57) 




(58) 



Note that all these quantum-statistical moments can be separated from each other 
by using their dependences on the phase shift tp (for the measurement of the corre- 
sponding spectral properties, see Vogel, W., [1995]). 
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2.2. Heterodyne detection 

It is worth noting that multiport homodyning for measurement of the complex ampli- 
tude is equivalent to (four-port) heterodyne detection (Yuen and Shapiro, J.H., [1980]; 
Yuen [1982]; Yuen and Chan [1983]; Shapiro, J.H., [1985]; Shapiro, J.H., and Wagner 
[1984]; Caves and Drummond [1994]). In the scheme, an optical field is combined, 
through a beam splitter, on the surface of a photodetector with a strong local oscilla- 
tor whose frequency Uq is offset by an amount Aoo from that of the signal mode in the 
optical field (Auj<^ujo). The measured photocurrent is electrically filtered in order 
to select the complex valued component at frequency Au. The classical statistics of 
this component correspond, under certain circumstances, to the quantum statistics 
of the two-mode operator 

a = as + a\, (59) 

where the subscripts S and /, respectively, are used to denote the signal mode at 
frequency ujq + Au; and the imaging mode at frequency c^o — Auo. Obviously, the 
imaging mode can be used to probe the signal mode. To show this, we first note that 
the measured quantity a can always be brought into the form 

a = &i + ia 2 , (60) 

where a\ and &2 are commuting Hermitian operators having a quantum-mechani- 
cal joint probability density p(a) =p(ai,a2), which can be found easily through its 
characteristic function 

(f)((3) = (e^-P*^ = J d 2 ae^-^>(a). (61) 

In particular, if the imaging mode is in the vacuum, then measurement of p(a) can 
be seen to yield the Q function of the signal mode, because of 

<!>(P) = (e-P &8 eP & s} = $(p-,s = -l) (62) 

(cf. Appendix B.4.). 

2.3. Parametric amplification 

Finally it should be mentioned that linear amplification of a signal mode and mea- 
surement of the complex amplitude of the amplified signal may be regarded as an 
equivalent of the homodyne measurement of the complex amplitude of the original 
signal (Bandilla and Paul, H., [1969, 1970]; Paul, H., [1974]). If the signal mode is 
strongly enhanced such that it behaves classically, the conjugated quadrature compo- 
nents can be measured simultaneously. The scheme can be realized by (nondegener- 
ate) parametric amplification, in which the photon destruction operators a' s and a'j, 
respectively, of the amplified signal and idler modes are obtained from the SU(1,1) 
input-output relations 

a's = Vg a s + y/g-l a\ (63) 



a 'i = \jd -la s + Vga\, 



(64) 
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where g is the (linear) signal-gain factor (for details, see, e.g., Mollow and Glauber 
[1967a,b]; Caves and Schumaker [1985]; Yurke, McCall and Klauder [1986]). 

It can be shown that when the signal mode is mixed with an idler mode that is 
initially in the vacuum state and the strong-amplification limit is realized, i.e., g — > 
oo, then the Q function of the original signal mode can be inferred from the measured 
distribution function for the complex amplitude of the amplified signal mode by ap- 
propriately rescaling this distribution (Leonhardt [1994]; see also Leonhardt and Paul, 
H., [1995]). It is worth noting that a degenerate parametric amplifier that operates 
as a squeezer [a\ — * e l ^a\ in eq. (63)] may be used to (partially) compensate for the 
detection losses in balanced four-port homodyne detection for measuring quadrature 
components (Leonhardt and Paul, H., [1994a]; § 3.9.1.), since it makes it possible to 
amplify a quadrature component with no increase of noise (Caves [1982]). 

Further, parametric amplification can be used to determine the quantum state of 
the signal mode by direct photodetection of the amplified signal if the idler mode is 
initially prepared in, e.g., a coherent state, without restriction to the Q function and 
without restriction to the strong-amplification limit (Kim, M.S., [1997a]; § 3.5.). In 
particular, when the idler mode is prepared in a strong coherent state such that it 
can be treated classically, ai —> ai, then - with regard to the photocounting formula 
(C.5) for detecting the amplified signal - eq. (63) takes the form of eq. (48), with 
a — —[(g — I)/ ' gY^a*! in place of eq. (49). In this case, eq. (50) also applies to the 
detection of the amplified signal, where now h = h s and r] = gr] D . 23 

2.4. Measurement of cavity fields 

Let us now consider the possibilities of phase-sensitive measuements of high-Q cav- 
ity fields. The realization of such fields, for example in a micromaser (Meschede, 
Walther and Miiller [1985]; Rempe, Walther and Klein [1987]; Brune, Raimond, Goy, 
Davidovich and Haroche [1987]), necessitates a strong suppression of cavity losses. 
That is, only a very small fraction of the field escapes from the cavity and could be 
used, e.g., in homodyne detection. Since the quantum efficiency of such a scheme is 
very low, it is impossible, in general, to extract the (complete) quantum statistics 
of the field inside the cavity from the homodyne data measured outside the cavity. 
To overcome the problem, methods have been developed which use atoms that travel 
through the cavity and probe the field inside. The system that is directly observed 
in the experiment is not the cavity field itself but the atoms after their interaction 
with the field. Under certain circumstances, the measured atomic occupation prob- 
abilities can then be related uniquely (in a more or less direct way) to the quantum 
state of the cavity field (§ 3.6.). Let us consider a measurement scheme typically 
used in Rydberg-atom superconducting cavity QED (see, e.g., Brune and Haroche 
[1994] and references therein). 24 The atoms that cross the cavity are prepared, after 

23 Note that eq. (63) implies that the signal and idler operators in the photocounting formula (C.5) 
are subject to different orderings, so that eq. (48) does not apply in general. 

24 Here circular Rydberg atoms are used in which one electron is placed on a highly excited energy 
level with large value of principal quantum number n, and the angular momentum takes its maximum 
value (I = m = n — 1). For n w 50 the atomic-dipole decay time is of the order of magnitude of 10~ 2 s. 
The cavity damping time is of the order of magnitude of 1CT 3 . . . 10 _1 s for superconducting cavities 
at subKelvin temperatures. The Q factor defined as the ratio of the cavity frequency to the cavity 
damping rate (inverse damping time) then reaches values of about 10 8 . . . 10 10 . The Rabi frequency 
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Figure 7: Set-up for the detection of a high-Q cavity mode. The Rydberg 
atoms are prepared in the state |e) in the box CB and cross a microwave 
zone R~ (source S2) between CB and C which mixes the (excited) state |e) 
and the (ground) state \g) before the atoms enter the cavity C. A second 
microwave pulse (source S2) applied between C and D in the zone R + again 
mixes the states |e) and \g) before the atoms reach the detector D, which 
renders it possible to detect linear superpositions of these two states. The 
microwave source Si can be used to feed into the cavity a (classical) field. 
(After Brune and Haroche [1994].) 



leaving an oven and before entering the cavity, into a superposition of Rydberg states 
|e) and \g) of different energies (Fig. 7). After preparation of the state |e), which 
involves laser and radiofrequency excitation, a microwave (Ramsey) zone can be used 
to apply a resonant pulse to the atoms, mixing the state |e) with another state \g) of 
different energy. In this way an initial atomic state of the form c e |e) + c g \g) can be 
injected into the cavity. The atoms then interact with the cavity mode while crossing 
it. After leaving the cavity, the atoms can again cross a microwave zone in order to 
realize state mixing before they enter a detector for measuring their energies by field 
ionization. This detection can be made energy-sensitive and one can thus count the 
atoms in |e) and \g). The dependence of the measured data on the various parame- 
ters characterizing the experimental conditions (such as state mixing and interaction 
time) can then be used to gain insight into the quantum state of the cavity mode 
(Vogel, W., Welsch and Leine [1987]; Bardroff, Mayr and Schleich [1995]; Bardroff, 
Mayr, Schleich, Domokos, Brune, Raimond and Haroche [1996]). 

The cavity and the atomic transition can be tuned to coincide (resonant regime) 
or slightly detuned (dispersive regime), with a frequency difference 5u> such that 
any exchange of energy between atom and field is made impossible. In the resonance 
regime, the system can be described (in its simplest version) by the Jaynes-Cummings 
model (Jaynes and Cummings [1963]; Paul, H., [1963]; for a review, see, e.g., Shore 
and Knight [1993]). The Hamiltonian consists of three terms describing the free two- 
level (circular Rydberg) atom, the free cavity mode and the atom-field coupling, the 
latter being given by 

H' = hn (a + a + a^a^j . (65) 

Here, a + = \e)(g\ and <r_ = |<?)(e|, respectively, are the atomic raising and lowering 
operators, and k is the coupling constant between the atom and the cavity mode 

which is typically of the order of magnitude of 10 . . . 100 kHz, is large compared with the atomic and 
the cavity decay rates, so that relaxation can be disregarded in first approximation (strong coupling 
regime) . 
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(Q = 2k, is the vacuum Rabi frequency, which is the rate at which the atom and the 
empty cavity exchange a photon at exact resonance.) 

In the dispersive regime, the atom-field interaction produces a dephasing of the 
field and also dephases the atom's wave function by an angle depending upon the 
number of photons in the cavity and on the quantum state of the atom. The nonres- 
onant interaction can be described by an effective Hamiltonian, 



The phase difference induced by the interaction of the two parts of the system can be 
measured when the two microwave zones in Fig. 7 are active. It can be shown that 
the probability that a probe atom which is initially prepared in the state |e) is finally 
detected in the state \g) (i.e., after the combined action of the two microwave zones on 
it) oscillates as a function of the microwave frequency, the position of the intereference 
(Ramsey) fringes being dependent on the additional dephasing produced by the cavity 
mode. It is worth noting that the Ramsey interferometer outlined here can serve as an 
apparatus for a quantum non- demolition measurement (QND) (Braginsky, Vorontsov 
and Khalili [1977]; also see Braginsky and Khalili [1992]) of the photon number in 
the cavity (Brune, Haroche, Lefevre, Raimond and Zagury [1990]; Brune, Haroche, 
Raimond, Davidovich and Zagury [1992]). 

When in the dispersive regime an additional microwave generator is used in order 
to resonantly couple a classical (strong) oscillator to the cavity mode, then the Wigner 
function of the cavity mode can, in principle, be measured directly (Lutterbach and 
Davidovich [1997]). Owing to the coupling of the classical oscillator to cavity mode, 
a displacement in phase space of the (initial) quantum state of the cavity mode is 
performed such that the density operator g is replaced with D^(a)gD(a). 25 The 
atoms initially prepared in state |e) cross the first Ramsey zone where they see a 7r/2 
pulse, so that |e) — > 2~ 1 / 2 (|e) + e' l ^\g)). The atoms then interact dispersively with 
the cavity field, and after that they undergoe a 7r/2 pulse in the second Ramsey zone 
before detection. Let P e and P g be the probabilities of detecting a probe atom in the 
states |e) and \g), respectively. It can be shown that the difference between these 
probabilities is 



where 0= (Q 2 /5)r is the additional dephasing of the atoms in the cavity, with r being 
the interaction time. In particular, if = ?/> = 7r/2, then 



which reveals that indeed the Wigner function is detected, provided that during the 
(ensemble) measurement the phases are controlled and each probe atom meets the 
cavity field in exactly the same state. 

The method of displaced states 26 can also be used for quantum-state measurement 
in the resonance regime, without phase- sensitive state mixing (Wilkens and Meystre 

25 Note that this coupling acts in a similar way as the mode mixing at an unbalanced beam splitter 
outlined in § 2.1.7.. 

26 An analogous method can be used in order to determine the quantum state of the center-of-mass 
motion of trapped ions (cf. § 4.2.2.). 




(66) 



AP = P e - P g = -Re [e l ^ + ^ (p{a) e 2i ^' a b\a))} , 



(67) 



AP = \nW{a), 



(68) 
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[1991]; Bodendorf, Antesberger, Kim, M.S., and Walther [1998]). Here, the probe 
atoms are prepared in one of the two considered states |e) and \g) or in a statistical 
mixture of them. The measured two-level occupation statistics of the atoms as a func- 
tion of the interaction time and/or the displacement parameter contain all knowable 
information on the quantum state of the cavity mode. 

The method of probing a cavity field using two-level atoms that travel through 
the cavity can also be extended to other than purely two-level probe atoms. So 
it was proposed to completely transfer, by adiabatic passage, the quantum state of 
the field onto an internal Zeeman submanifold of an atom (Walser, Cirac and Zoller 
[1996]). Utilizing a tomography of atomic angular momentum states by Stern-Gerlach 
measurements (Newton and Young [1968]), this angular momentum state can then 
be determined uniquely, by a finite number of magnetic dipole measurements. 

The spatial variation of the cavity mode implies a dependence on the atomic 
position of the Rabi frequency, which offers the possibility of probing the quantum 
state of the mode via atomic deflection (Akulin, Fam and Schleich [1991]; Herkommer, 
Akulin and Schleich [1992]; Freyberger and Herkommer [1994]). In the scheme in 
Fig. 8 a narrow beam of resonant atoms is proposed to pass through a node region of 
the mode, where the interaction Hamiltonian can be assumed to be a linear function 
of the atomic position. The cavity photons repulse the atoms so that the transverse 
momentum distribution of the atoms is changed during their passage through the 
cavity. When the cavity mode is (initially) in a pure state, then the measurable 
change provides (for appropriately prepared atoms) all the information needed for 
reconstructing it. 



class, field 
slit 




Figure 8: Atomic-deflection measurement scheme. A classical field prepares 
the two-level atoms in a superposition \g) + e llp \e) of the ground state \g) 
and the excited state |e). The so-prepared atoms pass a narrow slit which 
confines them to a region Ax <C A centred around the node of the standing 
light field at x = 0. After the resonant interaction with the field the deflected 
atoms travel to a screen put up far away from the resonator. The spatial 
distribution on this screen reflects the momentum distribution of the atoms 
right after the interaction. (After Freyberger and Herkommer [1994].) 



To extend the method to the reconstruction of mixed states, the effects of the 
classical field and the cavity field on the atoms can be combined such that the classical 
field travels orthogonally with respect to the cavity mode and both fields overlap in the 
region where the atoms cross them (Schneider, Herkommer, Leonhardt and Schleich 
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[1997]). In particular, when the classical field is sufficiently strong, then it can play 
a similar role as the local oscillator in balanced four-port optical homodyning. It is 
worth noting that under certain circumstances the transverse momentum distribution 
of the atoms which is observed in the dispersive regime is nothing but a smeared 
quadrature-component distribution of the cavity mode (§ 3.6.2.). 

3. Quantum-state reconstruction 

When in an experiment a sequence of quantities A\ {A\ = A\) can be measured, such 
that the density operator of the sytem under study can be given by 

£ = (69) 

i 

then all knowable information on the quantum state of the system can be obtained 
in principle (Fano [1957]). 27 In practice it must be ensured that enough data can 
be measured and that they can be decoded with reasonable accuracy and acceptable 
effort - a task which is not trivial in general. In particular if the measured data are 
directly related to a smeared quantum-state representation in the phase space, it may 
be an effort to derive from it the quantum state in another basis for studying spe- 
cific quantum-statistical properties. Typical examples of complete and overcomplete 
basis sets of observables for which measurement schemes have been designed are the 
generalized projectors M — > II(x, p>) = 7r _1 tp)(x, (p\ [see eq. (B.31)] and Ai — > 11(a) 
= n~ 1 \a)(a\ [see eq. (B.20) with s — — 1], respectively, for the quadrature-component 
states and the coherent states. If a scheme is not designed for measuring at least a 
[in the sense of eq. (69)] complete set of observables, then some a priori knowledge 
of the quantum state is required in order to compensate for the lack of data. 

3.1. Optical homodyne tomography 

From § 2.1.2. we know that in a succession of (phase-shifted) balanced four-port 
homodyne measurements the quadrature- component distribution p(x, ip) of a signal 
mode can be obtained for various values of the phase <p, provided that 100% quantum 
efficiency is realized. The quantum state is then known when p(x, ip) is known for 
all values of tp within a n interval (Vogel, K., and Risken [1989]). That is to say, 
all quantum-statistical properties can be obtained from the quadrature-component 
distributions measured in a ir interval. 

In particular, p(x, ip) can be used to reconstruct the Wigner function. From 
eq. (B.28) (s = 0) it can be derived that 

p(x : <p) — J dyW(xcosp — ysinp,xsinp+ycosp), (70) 

27 In connection with the study of the Pauli problem (Pauli [1933]) the sufficient set of observables 
that has to be mesurcd in order to obtain the complete information on the quantum state was also 
called quorum of a quantum-state measurement (Park and Band [1971]). It was shown that for the 
one-dimensional motion of a particle a quorum can by given by all symmetrized products of powers 
of position and momentum (cf. § 3.8.1.) or alternatively, all time derivatives of all powers of the 
position (Band and Park [1979]; Park, Band and Yourgrau [1980]). 



32 



3. QUANTUM-STATE RECONSTRUCTION 



where the definition W(q,p) =2~ 1 W[a = 2~ 1 ^ 2 (q + ip)] is used. Equation (70) reveals 
that the quadrature-component distributions can be regarded as marginals of the 
Wigner function. An integral relation of the form given in eq. (70) is also called 
Radon transformation (Radon [1917]). Inversion of the Radon transformation (70) 
yields the Wigner function in terms of the quadrature-component distribution. 28 From 
eq. (B.29) (s = 0) the Wigner function is derived to be 29 

W{q,p) = ^— ^ J dtp J dz J dx \z\ exp[iz(q cos tp+p sin tp — x)]p(x, tp). (71) 

Performing the z integral first would lead to an integral kernel that is not well- 
behaved. To overcome this difficulty in the numerical calculation, regularization 
techniques can be used. In the filtered back projection algorithm, the z integral is 
truncated such that 

W(q,p) ~ W c (q,p) = J dtp J dxK c (qcos tp+psin tp-x) p(x, tp), (72) 



K c {y) = \ dz\z\e iyz (73) 



with 



(z c > 0). Here, rapid variations of the quadrature-component distributions which 
correspond to frequencies higher than the cut-off frequency z c are suppressed, i.e., 
the exact distributions are effectively replaced with somewhat smeared ones. 

In the pioneering experimental demonstration of the method also called optical 
homodyne tomography (Smithey, Beck, Raymer and Faridani [1993]; Smithey, Beck, 
Cooper, Raymer and Faridani [1993]), a pulsed signal field is superimposed with a 
pulsed coherent-state field much stronger than the signal. The quadrature- component 
distributions are measured for a squeezed signal field and a vacuum signal field. The 
squeezed field is generated by using a walk-off compensated, travelling-wave optical 
parametric amplifier. The generated down-conversion signal centred at 1064 nm con- 
sists of two orthogonally polarized fields, the signal and idler, and has a bandwidth 
estimated to be 10 4 times that of the laser pump (532 nm, 300 ps, 420 pulses per 
second). The laser pump and the local-oscillator field (1064 nm, 400 ps) are obtained 
from a common laser source, and each local-oscillator pulse contains a mean number 
of photons of about 4 x 10 6 . The interfering fields are detected with high-quantum- 
efficiency (~ 85%) photodiodes, and the resulting current pulses are temporally in- 
tegrated and subtracted. The measurements and reconstructions are performed for 
a squeezed signal field and for a vacuum signal field, Figs. 9 and 10, the field mode 
detected being selected by the spatial-temporal mode of the local-oscillator field. The 
method of optical homodyne tomography was also extended subsequently to the con- 
tinuous wave-regime, including squeezed vacuum with a high degree of quantum noise 
reduction (Breitenbach, Muller, Pereira, Poizat, Schiller and Mlynek [1995]; Schiller, 
Breitenbach, Pereira, Muller and Mlynek [1996]) and bright squeezed light (Breiten- 
bach and Schiller [1997]; Breitenbach, Schiller and Mlynek [1997]). 



28 Inverse Radon transformation techniques are well known from tomographic imaging (for math- 
ematical details, see, e.g., Natterer [1986]). In quantum mechanics, the problem of tomographic 
reconstruction of the Wigner function of a particle moving in one dimension was first addressed by 
Bertrand and Bertrand [1987]. 

29 For phase-randomized states, see Leonhardt and Jex [1994]. 
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Figure 9: (a) In balanced four-port homodyne detection measured quadra- 
ture-component distributions at various values of the local-oscillator phase 
[P^Xy) corresponds to p(x,<p)]. (b) Variances of quadrature components 
vs local-oscillator phase: circles, squeezed state; triangle, vacuum state. In 
the experiment 4000 repeated measurements of the photoelectron difference 
number at 27 values of the relative phase tp are made. (After Smithey, Beck, 
Raymer and Faridani [1993].) 



As mentioned in § 2.1.2., the measured quadrature-component distributions do 
not correspond, in general, to the true signal mode, but they must be regarded as the 
distributions of a superposition of the signal and an additional noise source [eq. (21)], 
because of non-perfect detection. Substituting in eq. (70) for p(x, tp) the measured 
distributions p(x,tp;rj) with i]< 1 [eq. (19)] and performing the inverse Radon trans- 
form on them, the Wigner function of a noise-assisted signal field is effectively re- 
constructed. Equivalently, the reconstructed Wigner function can be regarded as an 
s-parametrized phase-space function of the true signal field, however with s<0. The 
characteristic function ^(z, tp; rj) of p(x, tp; rj) typically measured when 77 < 1 and the 
characteristic function s) of the phase-space function P(a; s) are related to each 
other according to eq. (B.27), with ty(z,tp;r)) in place of ty(z,tp) and s — 1 -h 77 1 in 
place of s in the exponential [cf. eq. (23) and footnote 7]. Hence, making in the 
exponentials in eqs. (B.28) and (B.29) for s the substitution s — 1 + r/~ 1 yields the re- 
lations between P(a; s) and p(x, tp; rf). In particular, eq. (B.29) can be used to obtain 
any phase-space function P(q,p; s) = 2 _1 P[a = 2 _1 / 2 (g + ip; s)] from the measured 
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quadrature-component distributions in principle: 30 

P(q,p;s) = — dp J dz J dx | exp ^(s-l+^^z 2 

x \z\ exp[iz(q cos ip+p simp — x)]p(x, (p; 77) j. (74) 

Obviously, when s= 1 — 7? -1 , then eq. (74) takes the form of eq. (71); i.e., replacing in 
eq. (70) p(x,ip) with p(x,ip;r)), inverse Radon transformation yields the signal-mode 
phase-space function P(q,p; s — 1 —i]' 1 ) in place of the Wigner function, 

p(x,(p;rf) — J dy P(x cosp—y sinip, x smtp+y cosp; s = l — rj^ 1 ). (75) 

The proposal was made to combine squeezing and balanced homodyning such that 
generalized quadrature-components 

x([/,,v) = fxq + up (76) 

for all real parameters \x and v can be measured, and to reconstruct the quantum state 
from the corresponding distributions p(X,/j,, v) (Mancini, Man'ko, V.I., and Tombesi 
[1995]; D'Ariano, Mancini, Man'ko, V.I., and Tombesi [1996]; Man'ko, O.V., [1996]; 
Mancini, Man'ko, V.I., and Tombesi [1997]). 31 In fact, x(/i, v) can be related to x(p) 
since it represents a scaled quadrature-component, 

x(fi, v) = \F\ (ae-^ + a)e itfi ) = 2 1 / 2 |F| %>), (77) 

with 

2 1 /2| F | = ^J^ + u 2 , tany? = - (78) 
(cf. § 2.1.2. and Appendix B.), which imples that 

p(x, n, v) = -j=±== p{-j^==^ = ^ ^ . (79) 

Hence, measurement of a particularly scaled quadrature-component distribution by 
means of an ordinary homodyne detector already yields all scaled quadrature com- 
ponents. Performing the analysis with a variable scaling parameter, it can be shown 
that eqs. (70) and (71), respectively, can be given in terms of p(x, /x, v) by 32 

p(x,/i,iy) = 7^f dk f d( lf dpe- lk{x - q »- pu) W{q,p) (80) 



30 When s > 1 — ry -1 , then in eq. (74) an inverse Gaussian occurs which may lead to an artifical 
enhancement of the inaccuracies of the measured data, so that a stable deconvolution might be 
impossible and the noise dominates the reconstruction of P(q,p; s). This effect is not observed when 
s < 1 — rj^ 1 , and hence a stable reconstruction with reasonable precision of P(a; s) for s < 1 — r/^ 1 
may therefore be expected to be feasible for any quantum state (see also § 3.9.). In particular, 
reconstruction of the Q function is always possible if 77 > 1/2. When s< 1 — i]^ 1 , then in eq. (74) 
the z integral can be performed first to obtain P(q,p; s) in a form suited for application of sampling 
techniques (§ 3.3.1.): P(q,p;s) = f£ dp J dx K(q,p 7 x,ip; s; 77) p(x , <p; rj) , with K(q,p, x, ip; s; 77) being 
a well-behaved integral kernel. 

31 For an application of the method (also called symplectic tomography) to trapped-ion quantum 
state reconstruction, see Mancini, Man'ko, V.I. and Tombesi [1996]; Man'ko, O.V. [1997]. 

32 For a detailed discussion of the transformation properties, see also Wiinsche [1997]. 
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Figure 10: From the measured quadrature-component distributions recon- 
structed Wigner distributions for (a),(b) a squeezed state and (c),(d) a 
vacuum state, viewed in 3D and as contour plots, with equal numbers of 
constant-height contours [W(X,P) corresponds to W{q,pj\. The recon- 
struction is performed by using inverse Radon transformation according to 
eq. (72). (After Smithey, Beck, Raymer and Faridani [1993].) 



and 

2 

W(q,p) = JdxjdfiJ dve iz{ - x -^- p ^p{x,n,v). (81) 

Comparing eq. (81) with eq. (71), we see that in any case a three-fold Fourier trans- 
formation is required in order to obtain the Wigner function from the homodyne data. 
Note that z in eq. (81) can be chosen arbitrarily, which reflects the above mentioned 
fact of overcomplete data. 



3.2. Density matrix in quadrature-component bases 

The reconstructed Wigner function can be used in order to calculate the density 
matrix in a quadrature-component basis (Smithey, Beck, Raymer and Faridani [1993]; 
Smithey, Beck, Cooper, Raymer and Faridani [1993]). The definition of the Wigner 
function given in Appendix B.4. can be rewritten, on expanding the density operator 
in a quadrature-component basis, as 

W(q,p) = - I dx'e 2iyx '(x-x',ip\g\x+x',ip), (82) 

7T J 



with 



x = q cos if + p sin </?, y = — q sin (p + p cos </?, 



(83) 
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which for tp — is nothing but the well-known Wigner formula (Wigner [1932]). 
Hence the density matrix in a quadrature- component basis can be obtained by Fourier 
transforming the Wigner function: 

(x — x', tp\g\x+x', tp) — J dy e~ 2tx y W(x cosy? — y sin tp, x sin tp+ycos tp). (84) 

Note that eq. (84) reduces to eq. (70) for x' = 0. 

It is worth noting that the reconstruction of the density matrix from the homodyne 
data can be accomplished with two Fourier integrals, avoiding the detour via the 
Wigner function (Kuhn, Welsch and Vogel, W., [1994]; Vogel, W., and Welsch [1994]). 
Writing the density-matrix elements as 

(x-x',tp\g\x+x',tp) = -!- / dze~ ixz ^(z,x', tp), (85) 

2lX J 

the characteristic function ty(z,x', tp) can be shown to be the characteristic function 
of a quadrature-component distribution, 

■*(z,x',tp) = y(z,tp), (86) 

with 

z = z(z,x') = z 2 + (2x') 2 and tp = tp - \it + axg(2x'+iz) . (87) 

Hence, the density-matrix elements can be obtained from the quadrature- component 
distributions by means of a simple two-fold Fourier transformation: 33 

(x-x',tp\g\x+x',tp) — -!- J dze~ txz J dxe tzx p(x, tp). 

It is worth noting that eq. (88) can be used to obtain the density matrix in different 
representations, by varying the phase of the quadrature component defining the basis. 
Further, eq. (88) can also be extended, in principle, to imperfect detection, expressing 
ty(z, tp) in eq. (85) in terms of ^(z, tp; rj) (cf. footnote 7; for details, see Kuhn, Welsch 
and Vogel, W., [1994]; Vogel, W., and Welsch [1994]). 

For the numerical implementation of the reconstruction based on eq. (88) spline- 
expansion techniques can be used (Zucchetti, Vogel, W., Tasche and Welsch [1996]), 34 

(x-x',tp\g\x+x',tp) ~ ^2K mn (x,x',tp)p(x m+1 ,<p n+1 ), (89) 

m,n 

where 

K mn (x,x',tp) = (27r)" 1 ^5 A£mim [5(^ +1 ,x , )] J B A2nin (^ +1 )5 ACfcifc (-a;), (90) 

k 

with z n = — 2x' cot(y5 n — tp) [Bf\(k), Fourier transform of B^(x)]. As an illustration 
of this method, in Figs. 11(a) and 11(b) the reconstructed density matrices in the 



33 For the two-fold Fourier transformation that relates the density-matrix elements to p(x,fj,,v), 
eq. (80), see D'Ariano, Mancini, Man'ko and Tombesi [1996]. 

34 Choosing a finite set of nodes {x n }, an approximate spline function f s (x) of f(x) is given by 
fs{x) = J2n f( x n+i)BAx„.n( x )i where B Axnin (x) = (x-x n )/Ax n if x n < x < x n+1 , B Axntn (x) = 
(x n +2 — x)/Ax n+ i if x n+ i <x< x„ +2 , and B Axn ,n{x) — elsewhere, and Ax n = x n+1 — x n ; for 
mathematical details of spline expansion, see de Boor [1987]. 
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"position" basis, <p = 0, and the "momentum" basis, (p = n/2, of a squeezed vacuum 
state are shown. The homodyne data were recorded by T. Coudreau, A.Z. Khoury 
and E. Giacobino. In the underlying experimental scheme the squeezing effect is 
obtained in a probe beam that interacts with cold atoms in a nearly single-ended 
cavity (Lambrecht, Coudreau, Steinberg and Giacobino [1996]). The phases tp — 
and ip = ir/2 that define the quadrature-component bases in Fig. 11 coincide with the 
phases of minimal and maximal field noise, respectively. 



Figure 11: From measured quadrature-component distributions recon- 
structed (real) density matrix (x, <f\g\x', <p) of a squeezed state in (a) the 
"position" basis (ip — 0) and (b) the "momentum" basis ((p — n/2). The ho- 
modyne data were obtained by T. Coudreau, A.Z. Khoury and E. Giacobino, 
using the experimental setup reported by Lambrecht, Coudreau, Steinberg 
and Giacobino [1996]. In the experiment, the quadrature- components were 
measured at 48 phases, and at each phase 7812 measurements were per- 
formed. 



3.3. Density matrix in the Fock basis 

Stimulated by the tomographic reconstruction of the Wigner function (§ 3.1.) much 
effort has been made to obtain the density matrix in the Fock basis from measurable 
data as direct as possible. Let us again start with analysing balanced four-port 
homodyning. 

3.3.1. Sampling of quadrature-components 

The problem of reconstruction of the density matrix in the Fock basis from the 
quadrature-component distributions can be solved, in principle, by relating the den- 
sity-matrix elements to derivatives of the Q function (cf. § 3.3.3.) and expressing 
the Q function in terms of the quadrature-component distributions, using eq. (B.29), 
with s — —l (D'Ariano, Macchiavello and Paris [1994a,b]). An equivalent formalism, 
which is suited for practice and which can also be applied to the reconstruction of the 
density matrix in other than the photon-number basis, is based on the expansion of 
the density operator as given in eq. (B.31). In the photon-number basis this equation 
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reads as 35 

Qmn = (m\g\n) = J dp J dx K mn (x, ip) p(x, ip). (91) 
The integral kernel (also called pattern function) 

K mn (x,<p) = {m\K(x,<p)\n) = f mn (x) S™-"* , (92) 

with 

K(x,ip) — — dz\z\ exp{iz [x(tp)— x]}, (93) 

2,71 J 

is studied in detail in a number of papers (D'Ariano [1995]; D'Ariano, Leonhardt 
and Paul, H., [1995]; Leonhardt, Paul, H., and D'Ariano [1995]; Leonhardt, Munroe, 
Kiss, Richter, Th., and Raymer [1996]; Richter, Th., [1996a]; Wunsche [1997]). The 
function f mn (x) (Fig. 12) is well behaved, and it is worth noting that it can be given 
by (Richter, Th., [1996a]; Leonhardt, Munroe, Kiss, Richter, Th., and Raymer [1996]) 

fmn(x) = ^- [lp m (x)(p n (x)] , (94) 

ax 

where ip m (x) and and <p m (x) are the regular and irregular (real) solutions of the 
harmonic-oscillator Schrodinger equation for a chosen energy value E m . 36 

Equation (91) reveals that the density matrix in the Fock basis can be sampled di- 
rectly from the measured quadrature-component statistics, since g mn can be regarded 
as a statistical average of the (bounded) sampling function K mn (x,ip) (D'Ariano, 
Leonhardt and Paul, H., [1995]; Leonhardt, Paul, H., and D'Ariano [1995]; Leon- 
hardt, Munroe, Kiss, Richter, Th., and Raymer [1996]; Leonhardt [1997c]; D'Ariano 
[1997a]). In an experiment each outcome x of x(<p), with ip e [0, 7r), contributes in- 
dividually to g mn , so that g mn is gradually building up during the data collection. 
That is to say, g mn can be sampled from a sufficiently large set of homodyne data in 
real time, and the mean value obtained from different experiments can be expected 
to be normal-Gaussian distributed around the true value, because of the central-limit 
theorem. Moreover, the sampling method can also be used to estimate the statistical 
error (see also § 3.9.1.). 

Experimentally, the method was successfully applied to the determination of the 
density matrix of squeezed light generated by a continuous- wave optical parametric 
amplifier (Schiller, Breitenbach, Pereira, Miiller and Mlynek [1996]; Breitenbach and 
Schiller [1997]; Breitenbach, Schiller and Mlynek [1997]). In the experiment the 
spectral component of the photocurrent in a small band around a radiofrequency Q 
is measured (overall quantum efficiency ~ 82%). In this case the measurement is 
on a two-mode quadrature-component x(ip) = e^ R [a{uo + Q)e~ llf + a(u> — Q)e lip ], <fin 
being the phase of the radio-frequency local oscillator, so that the scheme is basically 
a heterodyne detector. Examples of reconstructed diagonal density-matrix elements 
are shown in Fig. 13 (for reconstructed off-diagonal density-matrix elements, see, e.g., 
Schiller, Breitenbach, Pereira, Miiller and Mlynek [1996]; Breitenbach and Schiller 
[1997]). 



35 For the relation between the density-matrix elements and p(x,fj,,v), eq. (80), see D'Ariano, 
Mancini, Man'ko and Tombesi [1996]). 

36 Strictly speaking, the irregular (i.e., not normalizable) function ip n {x) has to be chosen such 
that ipnf'n ~ ip'nfn = 2/7r (for details, see Leonhardt [1997c]). Note that f mn {x) is not determined 
uniquely. 
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Figure 12: Plots of some typical pattern functions f mn (x) (lines) along with 
the products of regular wave functions ip n (x)ip m (x) (dashed lines). In (a) the 
diagonal pattern function with n = m = 4 is depicted. Oscillations between 
— 2/7r and +2/tt are clearly visible in the classically allowed region. Then 
the function swings over and decays like x~ 2 . In (b) the off-diagonal pattern 
function with n—1, m = 4: is depicted. It is less oscillating than the diagonal 
pattern function and it decays faster in the forbidden zone. In (c) and (d) 
highly oscillating pattern functions are shown for (c) n = m = 25 and (d) n 
= 20, m = 50. (After Leonhardt, Munroe, Kiss, Raymer and Richter, Th., 
[1996].) 



Further, the method was used successfully to measure the time-resolved photon- 
number statistics of a 5 ns pulsed field with a sampling time of 330 fs, set by the 
duration of the local-oscillator pulse (Munroe, Boggavarapu, Anderson and Raymer 
[1995]). From eqs. (91) and (92) it is easily seen that the photon-number probability 
distribution p n = g nn can be given by 

p n = 7T J dxf nn (x)p(x), (95) 

where p(x) = (27r) _1 / dipp(x, ip) is the phase-averaged quadrature-component distri- 
bution. Equation (95) reveals that for sampling the photon-number statistics the 
phase need not be controlled - a situation that is typically realized when the signal 
and the local oscillator come from different sources. In the experiment, an argon- 
laser-pumped Ti:saphire laser is used in combination with a chirped-pulse regenera- 
tive amplifier to generate ultrashort, transform limited local-oscillator pulses (330 fs) 
at a wavelength of 830 nm and a repitition rate of 4 kHz with approximately 10 6 
photons per pulse. The signal is from a laser diode of wavelength 830 nm and pulse 
width 5 ns. Figure 14 shows examples of the (with ~ 65% overall quantum effi- 
ciency) measured quadrature-component probability distributions and the resulting 
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Figure 13: Photon- number distribution of a squeezed vacuum and the vac- 
uum state (inset) reconstructed from the quadrature-component distributions 
according to eq. (91). Solid points refer to experimental data, histograms 
to theory. The experimentally determined statistical error is 0.03. (After 
Schiller, Breitenbach, Pereira, Midler and Mlynek [1996].) 



photon-number probability distributions at two times in the signal pulse. In particu- 
lar, the photon-number statistics are seen to change from nearly Poissonian statistics 
(laser above threshold) to thermal-like statistics (laser below threshold) (for ultrafast 
homodyne detection of two-time photon-number correlations, see § 3.8.1.). 

So far in the formulas, perfect detection has been assumed. The problem of ex- 
tending eq. (91) to imperfect detection such that p(x, ip) and K mn (x, if), respectively, 
are replaced with p(x,ip;r)) and a sampling function K mn (x,ip;ri) that compensates 
for the losses has also been considered, 37 

K mn (x, ip;r]) = (m\K(x, <p;rj)\n), (96) 

where 

K(x,<p;rj) = — J dz\z\ exp{iz [x(tp)-x\ + \{ri' 1 - l)^ 2 } (97) 

(D'Ariano [1995]; D'Ariano, Leonhardt and Paul, H., [1995]; Leonhardt, Paul, H., 
and D'Ariano [1995]; D'Ariano [1997a]; D'Ariano and Paris [1997a]). It has been 
shown that K mn (x, ip; rf) is a well-behaved bounded function provided that 77 > 1/2. 

It is worth noting that the reconstruction formula (91) also applies to other 
than harmonic-oscillator systems (Leonhardt and Raymer [1996]; Richter, Th., and 
Wunsche [1996a,b]; Krahmer and Leonhardt [1997b, 1997c]; Leonhardt [1997a]; Leon- 
hardt and Schneider [1997]). To be more specific, p(x,ip) in terms of g mn reads as 38 

p(x, v?) = E g mn (x)e-« m - n ^ g mn , (98) 

m,n 

37 Equation (96) follows from eqs. (B.31) and (B.32), expressing *$>(z, ip) in terms of ^(z, ip; 77); see 
footnote 7. 

38 Note that substituting in eq. (91) for p(x,ip) the expression on the right-hand side of eq. (98), 
carrying out the ip integral and using the orthonormalization relation (100) just yields an identity. 
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Figure 14: Measured phase-averaged quadrature-component distributions for 
(a) t = 4.0 ns and (b) t = 6.0 ns, and the the resulting time-resolved photon- 
number distributions for (c) £ = 4.0 ns and (d) t — 6.0 ns, obtained from (a) 
and (b) by using eq. (95). and p(n), respectively, correspond to p(x) 

and p n . (After Munroe, Boggavarapu, Anderson and Raymer [1995].) 



where 

9mn{x) = lp m (x) lj) n (x), (99) 

and from eqs. (91) and (98) together with eqs. (92), (94) and (99) it follows that 
the functions 7rf mn (x), eq. (94), are orthonormal to products of energy eigenfunctions 
g mn {x), eq. (99), 

n J dx f mn (x) g m ' n >(x) = 6 mm >6 nn > for E m -E n = E m , -E n > . (100) 

It can be shown that ip m (x) and (p m (x) can be the regular and irregular solutions, 
respectively, which solve a Schrodinger equation 

4>m{x) = E m (f) m (x) (101) 

for chosen energy E m , with U(x) being an arbitrary potential. 39 This offers the pos- 
sibility of reconstruction of the density matrix (in the energy representation) of a 
particle in an arbitrary one-dimensional potential U (x) from the time-dependent po- 
sition distribution of the particle. The quadrature-component distribution in eq. (91) 
has to be regarded as a time-dependent position distribution p(x,t) = (x\g(t)\x) and 
the phase integral converts into a time integral, 

7r r +T / 2 r 

Qmn = lim - / dt dx e tUmnt f mn {x) p(x, t), (102) 

[Qmn — Qmn(t)\t=o] v mni transition frequencies]. Obviously, the position distribution 
p(x, t) can also be used for reconstruction of other quantum-state representations, 



For a proof, see, e.g., Leonhardt and Schneider [1997]; Leonhardt [1997c]. 
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such as tomographic reconstruction of the Wigner function. Note that the time 
interval T, in which the position distribution can be measured is limited in general 
(for a data analysis scheme for determining the quantum state of a freely evolving 
one-dimensional wave packet, see Raymer [1997]; see also § 4.4.1.). 

In the numerical implementation of eq. (91) [or eq. (102)] the sampling function 
can be calculated on the basis of the analytical result given in eqs. (92) and (94), 
using appropriate numerical routines (Leonhardt, Munroe, Kiss, Richter, Th., and 
Raymer [1996]). An alternative way is the direct inversion of the underlying basic 
equation (98) that expresses the quadrature-component distributions p(x, tp) in terms 
of the density-matrix elements g mn (Tan [1997]). 40 For any physical quantum state 
the density-matrix elements g mn must eventually decrease to zero with increasing 
m(n). Therefore it follows that the expression on the right-hand side of eq. (98) 
can always be approximated to any desired degree of accuracy by setting g mn ps 
for m(n) >n max . To obtain the (finite) number of density-matrix elements from the 
measured quadrature-component distributions, the resulting equation can always be 
inverted numerically, using standard methods, such as least-squares inversion (Ap- 
pendix D.; see also § 3.9.2.). It should be pointed out that the method also applies 
when the detection efficiency is less than unity. Recalling eqs. (19) and (20), it is 
easily seen that p{x,ip;rj) is related to g mn according to eq. (98), with g mn (x) being 
replaced with 

9mn(x;r]) = J &x' i> m (x')il) n (x)p(x-x;vi). (103) 

When the characteristic function ty(z,ip) of the quadrature-component distribu- 
tion p(x,ip) can be measured directly (§ 4.2.1.), then the density matrix in the Fock 
basis can be obtained from ^(z,tp) by replacing in eqs. (91) and (92) the x integral 
with the z integral over f (z)^(z, </?), where 

l mn (z) = Jdxe^f mn (x). (104) 
More explicitly, the result can be given by (Wallentowitz and Vogel, W., [1996b]) 

r , -ik» /*°° i o(k), \ f Re ^( z ^f) if k even, 
Qnn +k = d(pe lk * dzSi k \z)\ 105) 
Jo Jo lm^(z, ip) if k odd 

(Jfc>0), with 



s? , (2)= ^/^:(£)' +i L (»» ( , 72)e -,vJ(- 2 ); /2 "-n, 

n y ' 7T \{n + k)\ \2j n v 1 ' \ (_2)( fc " 1 )/ 2 if k odd 

[L^\x), Laguerre polynomial]. 

So far, reconstruction of arbitrary quantum states has been considered, which can 
require measurement of the quadrature components at a large number of phases (see 
also §3.9.1.), which reveals that the Pauli problem [i.e., reconstructing a quantum 



40 The underlying basic equation for eq. (102) is given by eq. (98) with p(x, t) and e ~ lVmnt in place 
of p(x, if) and e ~ l ( m ~ n )v^ respectively. Inverting it numerically, the time need not be infinitely large 
as it might be suggested from the analytical result given in eq. (102) (Opatrny, Welsch and Vogel, 
W., [1997c]). 
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state from p(x, ip) and p(x, <p + 7r/2)] cannot be uniquely solved in general. However, 
if there is some a priori information on the quantum state to be reconstructed, then 
p(x, p) need not be known for all phases within a n interval. In particular when the 
state is known to be a pure state that is a finite superposition of Fock states, 

l*> = E c n\n), (107) 

n=0 

then it can be reconstructed from two quadrature-component distributions p(x, tp) and 
p(x, p + n/2). 41 In this case the problem reduces to solving blocks of linear equations 
for the unknown coefficients in the Fock-state expansion of the state (Orlowski and 
Paul, H. [1994]). 42 



3.3.2. Sampling of the displaced Fock-states on a circle 

From § 2.1.7. we know that in unbalanced homodyning the photon-number distri- 
bution of the transmitted signal mode is, under certain conditions, the displaced 
photon-number distribution of the signal mode, p m (a), the displacement parameter 
a = \a\e t,f being controlled by the local-oscillator complex amplitude [see eq. (51)]. 
Expanding the density operator in the Fock basis, p m {a) can be related to the den- 
sity matrix of the signal mode as 

Pm(oi) = (m,a\g\m,a) = ^2(n\m,a)(k\m,a)*g kn , (108) 

k,n 

where the expansion coefficients (n\m,a) can be taken from eq. (B.15) or eq. (B.16). 
Equation (108) can always be inverted in order to obtain g mn in terms of Pfc(a) 
(Mancini, Tombesi and Man'ko, V.I. [1997]; Mancini, Man'ko, V.I., and Tombesi 
[1997] ): 43 




(109) 



where 

K k mn {a) = (m\K k (a)\n) = (£±1)* (m\6(a-a; -s)\n) (110) 



41 The Pauli problem of determining the quantum state of a particle from the position distribution 
and the momentum distribution has been studied widely, and it has turned out that it cannot be 
solved uniquely even if the particle moves in a one-dimensional potential and is prepared in a pure 
but arbitrary state (for the problem including finite dimensional spin systems, see Pauli [1933]; 
Feenberg [1933]; Kemble [1937]; Reichenbach [1946]; Gale, Guth and Trammell [1968]; Lamb [1969]; 
Trammell [1969]; Band and Park [1970, 1971]; Park and Band [1971]; d'Espagnat [1976]; Kreinovitch 
[1977]; Corbett and Hurst [1978]; Prugovecki [1977]; Corbett and Hurst [1978]; Vogt [1978]; Band 
and Park [1979]; Park, Band and Yourgrau [1980]; Ivanovic [1981, 1983] Moroz [1983, 1984]; Royer 
[1985, 1989]; Friedman [1987]; Pavicic [1987]; Wiesbrock [1987]; Busch and Lahti [1989]; Woottcrs 
and Fields [1989]; Stulpc and Singer [1990]; Bohn [1991]; Weigert [1991]). 

42 These coefficients can also be obtained from p(x,ip) and d/d(p p(x, t£)\ v =o (Richtcr, Th., 
[1996c]). More generally, it can be shown that a pure state can always be determined from 
p(x,ip) and d / dip p(x , ip) , i.e., from the position distribution and its time derivative in the case 
of a particle that moves in a one-dimensional potenial and is prepared in an arbitrary pure state 
(Feenberg [1933], for the problem, also see Gale, Guth and Trammell [1968]; Royer [1989]). 

43 Equation (109) directly follows from eqs. (B.20) and (129). 
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is bounded for s E (—1,0], the 5 operator being given by eq. (B.19). This offers 
the possibility of direct sampling of g mn from the displaced Fock-state probability 
distribution p m (a). 

Since p m (a) as a function of a for chosen m already determines the quantum 
state (§§ 2.1.5. and 3.3.3.), it is clear that when m is allowed to be varying, then - 
in contrast to eq. (109) - p m {a) need not be known for all complex values of a in 
order to reconstruct the density- matrix elements g^ n from p m {a). In particular, it is 
sufficient to know p m (a) for all values of m and all phases if, \a\ being fixed (Leibfried, 
Meekhof, King, Monroe, Itano and Wineland [1996] 44 ; Opatrny and Welsch [1997]; 
Opatrny, Welsch, Wallentowitz and Vogel, W., [1997]). For chosen \a\ we regard 
p m (a) as a function of <p and introduce the Fourier coefficients 

Pm(l«l) = ^- F&pe^pja) (111) 

(s = 0,1,2,...), which are related to the density-matrix elements whose row and 
column indices differ by s as 

oo 

p s m (\a\) = Y,G S mn M)Qn + sn, (H2) 

n=0 

where 



G s mn M)=e-^m^n\(n+s)\ 

J L f_lV+IU|2(m+n-3-I)+s 

P ^ jKrn-mn+s-j)\U(m-l)\(n-l)\ ' 1 } 

with J = min(m, n + s) and L = min(m, n). Inverting eq. (112) for each value of s 
yields the sought density matrix. Since there has not been an analytical solution, 
eq. (112) has been inverted numerically, 45 setting f? m „~0 for m{n) > n max (cf. the 
last paragraph but two of § 3.3.1.) and using least-squares inversion (Appendix D.; 
§ 3.9.2.). In this way g m +sm can be given by 

N 

e m+sm = Y,Kn(W)p s n (\a\), (114) 

n=0 

where it is assumed that p n is measured for n — 0, 1, 2, . . . , N, with N < n max , and the 
matrix F^ n (\a\) is calculated numerically. Combining eqs. (114) and (111), Q m +sm 
can be given in a form suitable for statistical sampling. 

An extension of eq. (112) to nonperfect detection is straightforward. In this case 
eq. (50) applies, and the measured probability distribution P m (a; rj) can be related to 
p m (a) as shown in Appendix C From eqs. (Ill) and (C.3) it can be seen that when 
in eq. (112) the Fourier component is replaced with the actually measured 

one, then the matrix G^ n (|Q;|) has to be replaced with the matrix 

G s m Mv) = Y, P Mv)G s ln (\a\). (115) 
i 



44 Here the method was first used for reconstructing experimentally the density matrix of the 
center-of-mass motion of a trapped ion (also see § 4.2.2.). 
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For examples, see, e.g., Opatrny, Welsch, Wallentowitz and Vogel, W., [1997]. 
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Similarly, multichannel detection of the photon-number distribution can be taken into 
account. In particular when a photon-chopping scheme as outlined in § 2.1.7. is used, 
then the Fourier components of the measured iV-coincidence-event probability distri- 
bution can again be related to the density-matrix elements according to eq. (112), 
but now with 

G s mn (\a\, V , JV) = £ P mlk (N)P kll (r]) G s ln (\a\) (116) 
k,i 

in place of G^ n (|a|, 77), where P m \k(N) can be taken from eq. (53). 
3.3.3. Reconstruction from propensities 

As already mentioned, the displaced- Fock-state probability p m (a), eq. (108), as a 
function of a defines for each m a propensity prob(a) = n~ 1 p m ( y a) of the type given in 
eq. (40). The simplest example is the Q function (m = 0), which can be measured in 
perfect eight-port balanced homodyne detection, with the vacuum as quantum filter 
(§ 2.1.4.). Propensities contain all knowable information on the quantum state, and 
therefore all relevant properties of it can be obtained from them in principle. However, 
in order to obtain the "unfiltered" quantum state, the additional noise introduced by 
the filter must be "removed" , which may be an effort in practise. 

When the Q function is known, then the density matrix in the photon-number 
basis can be calculated straightfowardly, using the well-known relation 46 

n Qm+n 

Qmn = (m\g\n) = ===== - e H Q(a) 

\Jm\n\ oa* m oa n 

which corresponds to 

(118) 

a=a*=0 

(Wunsche [1991, 1996a]). Equation (118) can also be extended to other than vacuum 
filters, in order to obtain the density operator in terms of derivatives of more general 
propensities of the type given in eq. (40) (Wunsche and Buzek [1997]). In order to 
obtain the density matrix from the homodyne data measured, derivatives on them 
must be carried out, which probably could be done with sufficient accuracy only for 
states that contain very few photons. 

In practice it may be more convenient to handle integrals rather than derivatives. 
It can be shown that 47 




(119) 



where m>n and s < —1, L™(z) being the Laguerre polynomial (Paris [1996b]). In 
particular when s — 1 — 2n~ 1 , then P(a; s) is just the smoothed Q function measured 

46 Writing tt e^ 2 Q(a) — e'"' 2 (a|£|a) = J2 m n e'"' 2 (a\m) (m\g\n) (n\a) and recalling that (n\a) — 
(n!) -1 / 2 a n e~l Q l 2 / 2 , eq. (117) can be derived easily. 

47 Equation (119) can be derived, applying eq. (B.22) to F= \n)(m\ and calculating the c-number 
function F(a; s) associated with F in chosen order. Note that for s = —1 the integral form (119) 
corresponds to the differential form (117). 



(H7) 

'=rv*=n 



q — 7r exp 



a 
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in nonperfect detection (cf. § 2.1.4.). Unfortunately, eq. (119) is not suitable for 
statistical sampling, since the integral must be performed first and after that the 
summation can be carried out. Moreover, the inaccuracies of the measured P(a; s) 
together with the Laguerre polynomials can give rise to an error explosion in the 
reconstructed density matrix, so that an exact reconstruction of the density matrix 
from measured data may be expected to be possible only for states that contain finite 
(and not too large) numbers of photons. In this case the p sum in eq. (119) can be 
truncated at p = N — m, where the value N has to be chosen large enough to ensure 
that (a^dP) = for i,j > N. Now the p sum can be performed first, and a (state- 
dependent) integral kernel for statistical sampling can be calculated (Paris [1996b, 
1996c]). 

3.4. Multimode density matrices 

The extension of the methods outlined in §§ 3.2. and 3.3. to the reconstruction of 
multimode density matrices from the corresponding multimode joint quadrature- 
component distributions or multimode joint propensities is straightforward. The 
situation is somewhat different when combined distributions, i.e., distributions that 
are related to linear combinations of the modes, are measured. Let us consider the 
two-mode detection schemes shown in Fig. 4 in § 2.1.3.. When the sum quadrature- 
component distribution of two modes, ps(x, a, ifi, ip 2 ), is known for all phases ipi and 
<f 2 within 7r intervalls and all superposition parameters a G (0,7r/2), then it can be 
shown, on recalling eq. (30), that the reconstruction of the two-mode density matrix 
in a quadrature-component basis can be accomplished with a three-fold Fourier inte- 
gral (Opatrny, Welsch and Vogel, W., [1996, 1997b]): 

(xi — x' 1 ,X2 — x' 2 , (f\g\x 1 +x' 1 ,x 2 +x 2 , <p) 

= Q^j J d Zl J dz 2 e-^ +x ^ J dxe^ps(x,a,^^ 2 ), (120) 

where 

y = v /^ 2 + ^ + (2xi) 2 + (24) 2 , tana = 

and 

&k = <p - §7r + axg(2x' k +iz k ) (122) 

(A; = l,2). 48 

The generalization to the iV-mode case is straightforward. Suppose that we can 
measure the probability distribution of a weighted sum of quadratures x = X^Li 
fk({ai})xk, where ai are N — 1 parameters that can be controlled in the experiments, 
I — 1, 2, . . . N— 1, and fk are real functions which satisfy J2k=i fk({ a i}) = 1 identically 
for each set of parameters {a{\. Let gi be the inversions of fk, gi[{fk({®m})}] = 
ai (i.e., the set {a{\ parameterizes the surface of an iV-dimensional sphere). From 
the measured sum quadrature-component probability distribution ps(x, {a{\, {(fik}) 
its characteristic function ^ s (z, {ai}, {(fk}) can be calculated as a Fourier transform. 



-(24: 



z 2 1 + (2x' i y 



(121) 



For an extension to imperfect detection, see Opatrny, Welsch and Vogel, W., [1997b]. 
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The characteristic function of the joint quadrature-component probability distribution 
can then be calculated as 

*j({z k }, W k }) = ^ s (z, { 9l ({z k /z}, {<p k })}) , (123) 

[z = (J2k=i z "k) l ^ 2 ]i fr° m which the iV-mode density matrix can be obtained by an 
iV-fold Fourier transform; the whole iV-mode density matrix reconstruction is thus 
accomplished by an N+ 1-fold integration of the measured data. 

In the Fock basis the reconstruction of a two-mode density matrix from the com- 
bined quadrature-component distribution ps{x, a, ip±, </? 2 ) can be accomplished with 
a four-fold integration (Raymer, McAlister and Leonhardt [1996]; McAlister and 
Raymer [1997b]; Richter, Th., [1997a]): 

/f"K P7Y f'TT 

dx J^aJ^ dcp^ d</? 2 i?™(x,o;,</?i</?2)p s (x,Q;,</7i,</?2), (124) 

where the integral kernel i?™^^(x, a, </?i</? 2 ) is suitable for application of statistical 
sampling. It can be given by 

iCEfo a, ^ = C^(x,a)e i ^- Bl ^e i ^- na ^, (125) 

where 

CS; (. a) ^/d,/ mW ^^^),„(^). m 

Here, ip mi (y) and ip ni (y) are harmonic-oscillator energy eigenfunctions, and f m2 n 2 iu) 
is given by eq. (94). An alternative integral expression for reads 



C 2 n?(^a) = 77T^ ldz\z\e- lzx f (zcosa)f (zsina), (127) 

m\n\\ i ) (271") J ±-mini y ^£_m 2 n 2 v n v ' 

with / {£) being given by eq. (104). Since / (z) is the kernel function for recon- 
struction of the single-mode density matrix in the Fock basis from the quadrature- 
component characteristic function </>), from eqs. (105) and (106) it is seen that 
it can be expressed in terms of the associated Laguerre polynomial. Using this in 
eq. (127), then the integral can be performed to obtain a representation of r^^(x, a) 
as a finite sum over confluent hypergeometric functions (Richter, Th., [1997a]). 

Finally, the problem of reconstruction of the quantum state of unpolarized light 
was studied, by considering the (two-mode) density operator 

■k)(k,n-k^j (128) 

(Lehner, Leonhardt and Paul, H., [1996]). Obviously, the probability distribution p n 
of finding n completely unpolarized photons in the signal field characterizes uniquely 
the quantum state (128). Two schemes for determining p n were studied. First, ex- 
pressions for the sampling function were derived for the case when the two orthogonal 
polarization modes are delivered to two balanced homodyne detectors and the joint 
quadrature-component distributions are measured (Krahmer and Leonhardt [1997a]. 
Second, it was shown that measurement of the quadrature-component distributions 
of any linearly polarized component of the signal field is sufficient to determine p n , 
the corresponding sampling function being closely related to the single-mode function 
(Richter, Th., [1997b]). 
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3.5. Local reconstruction of P(a; s) 

The displaced photon-number statistics that is measurable in unbalanced homodyning 
(§ 2.1.7.) can be used for a pointwise reconstruction of s-parametrized phase-space 
functions (Wallentowitz and Vogel, W., [1996a]; Banaszek and Wodkiewicz [1996]). 49 
From eq. (B.21) together with eqs. (B.19) and (B.13) it is easily seen (cf. also Moya- 
Cessa and Knight [1993]) 

9 00 / <?4-1 \ 171 

p(q ' ; - 5) ^^,S(^t) p - (a) - (129) 

Hence, all the phase-space functions P(a; s), with s< 1, can be obtained fromp m (a;) 
for each phase-space point a in a very direct way, without integral transformations. 
In particular when s— —1, then eq. (129) reduces to the well-known result that Q(a) 
= 7r~ 1 po(tt), with po(a) = (a\g\a). Further, choosing s = in eq. (129), we arrive at 
the Wigner function, 

W» = -£(-l)> m (a). (130) 



n m=0 



Equation (130) reflects nothing but the well-known fact that the Wigner function is 
proportional to the expectation value of the displaced parity operator (Royer [1977]). 50 
Experimentally, the method was first applied to the reconstruction of the Wigner 
function of the center-of-mass motion of a trapped ion (Leibfried, Meekhof, King, 
Monroe, Itano and Wineland [1996]; § 4.2.2.). 

From § 2.1.7. we know that in unbalanced homodyning the quantum efficiency 
V = \Uki\ 2 VD is always less than unity even if r/o = 1, because of \Uki\ < 1. Equation 
(129) can be extended to imperfect detection in order to obtain (Wallentowitz and 
Vogel, W., [1996a]; Banaszek and Wodkiewicz [1997b]) 51 



77(s-l) + 2' 
rj(s-l) 



P m (131) 



[a = —(Uk2/Uki)aL, eq. (49)], where the measured photon-number distribution P m 
can be regarded as a smoothed displaced photon-number distribution p m (a; if) of the 
signal mode. The method works very well and is well suited for statistical sampling 
if s< l — rj^ 1 ; i.e., when the weighting factors improve the convergence of the series. 52 
Hence reconstruction of the Q function with reasonable precision is always possible if 
r]> 1/2 (for computer simulations of measurements, see Wallentowitz and Vogel, W., 
[1996a]; Banaszek and Wodkiewicz [1997b]). 

As already mentioned in § 2.3., measurement of the photon-number distribution 
of a linearly amplified signal can also be used - similarly to measurement of the 



49 For a squeezed coherent local oscillator, see Banaszek and Wodkiewicz [1998]. 
50 For proposals of measuring the Wigner function of a particle using this fact, see (Royer [1985, 
1989]). 

51 This can be easily proved correct, applying eq. (C.4) and expressing p m (a) in eq. (129) in terms 
of P m , eq. (50). 

52 The feasibility of reconstruction of the Wigner function of truncated states was also demon- 
strated for r/ < 1 (Wallentowitz and Vogel, W., [1996a]). 
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displaced photon-number statistics in unbalanced homodyning - for reconstructing 
the quantum state of the signal mode. It can be shown that when the idler mode is 
prepared in a coherent state, then the phase-space funcion P(a; s) of the signal mode 
can be related to the measured photon-number distribution as 53 



P(a;s) = 



2g 



n(2-g- 



E 

m=0 



r] D (s+g-2) + 2 
r) D (s+g-2) 



(132) 



(Kim, M.S., [1997a,b]), where a = -[(g-l)/g] 



1/2 



a, 



Note that the gain factor g and 



the quantum efficiency r] D of the detector enter separately into eq. (132) [in contrast 
to ?7= \U k i\ 2 n D in eq. (131)]. 



3.6. Reconstruction from test atoms in cavity QED 

Let us now turn to the problem of reconstruction of the quantum state of a high-Q 
cavity field from measurable properties of test atoms in detection schemes outlined in 
§ 2.4.. Though at a first glance the schemes look quite different from the homodyne 
detection schemes, there are a number of remarkable similarities between them. 



3.6.1. Quantum state endoscopy and related methods 

Let us first consider a two-level (test) atom that resonantly interacts with a single- 
mode cavity field according to a A;-photon Jaynes-Cummings model, the atom-field 
interaction Hamiltonian being given by 

H( k ) = (a + a k + a tfc <r_) , (133) 

which for k = 1 reduces to eq. (65). When the atoms are initially prepared in superpo- 
sition states |±) = 2 -1 / 2 (\g) ± |e)) and the excited-state occupation probabilities 
P e ± (t) are measured as functions of time, then the cavity- mode density-matrix ele- 
ments Qnn+k can be determined (Vogel, W., Welsch and Leine [1987]). To be more 
specific, it can be shown that the difference P^(t) — P e + (t) reads as 

oo 

P~(t) - P e + (t) = 2 ]T <#> sm(Q nn+k t) , (134) 

n=0 

where VL nn+k = 2k W {(n+1) ■ ■ ■ (n + k)} 1 / 2 , and 

a ( n ] = — Qnn+k + ex.; (135) 
li 

i.e., the off-diagonal density-matrix elements g nn +k can be obtained directly from 
the coefficients for two phases ip, such as ip = and ip — n/2. Provided that 
the interaction time t can be varied in a sufficiently large interval (0,T), the Fourier 

53 From eq. (C.5) it can be found that for a — Q, eq. (129) relates the measured distribution P m to 
the phase-space function of the detected mode at the origin of the phase space, Pdct(0, s). Equation 
(132) can then be proved correct, using the unitary tranformation (63) and expressing Pd c t(a, s) 
in terms of phase-space functions of the signal and idler modes by convolution (see also Leonhardt 
[1994]; Kim, M.S., and Imoto [1995]). 
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transform of P e (t) — P^(t) consits of sharp peaks, whose values yield the sought 
coefficients as 54 



a {k) 



= ^£dt sm(Q nn+k t) [P e -(t) - P e + (t)] (136) 

(T — > oo). To measure the diagonal density-matrix elements g nn , it is sufficient to 
prepare the atom in the excited state, P e (t)\ t= o = 1, and observe (for arbitrary k) the 
atomic-state inversion AP = P e — P g = 2P e — 1, 

oo 

AP(t) = Qnn cos(tt nn+k t) , (137) 

n=0 

from which g nn can be obtained by Fourier transformation. 55 

In cavity QED the 1-photon Jaynes-Cummings model is typically realized, so that 
the method - also called quantum state endoscopy - does not yield the off-diagonal 
density- matrix elements g nn +k with k > 1. When the quantum state is a priori known 
to be a pure state such that g mn is given by g mn = c m c* n with c m c* m+l ^OVm, 56 then 
eq. (134) [together with eq. (135)] can be taken at a sufficiently large number of time 
points (and at least at two phases) in order to obtain [after truncating the state at a 
sufficiently large photon number n max according to eq. (107)] a system of conditional 
equations for the expansion coefficients c m , which can be solved numerically (Bardroff, 
Mayr and Schleich [1995]; Bardroff, Mayr, Schleich, Domokos, Brune, Raimond and 
Haroche [1996]). 

The reconstruction problem for arbitrary quantum states can be solved by per- 
forming a displacement of the initial state of the cavity field such that g is replaced 
with (a)gD(a), and hence (k = l) 

oo 

AP(t) = Qnn(a) cos(pVn+lt) (138) 

n=0 

(Q = 2k, k = k^), where g nn (a) = (n\D\a)gD{a)\n) = (n,a\g\n,a). Now g nn {a) can 
again be obtained from AP(t) by Fourier transformation, and from g nn {ct) the quan- 
tum state of the cavity mode can be obtained, applying, e.g., the methods outlined in 
§§ 3.3.2. and 3.5.. Alternatively, the quantum state can also be reconstructed when 
the interaction time is left fixed and only a = \a\e t(p is varied (Bodendorf, Antesberger, 
Kim, M.S., and Walther [1998]). Regarding the measured atomic occupation proba- 
bilities as functions of a = \a\e lip , AP(t) = AP(t, a), and introducing for chosen t and 
| a | the Fourier coefficients 

1 r 2v 

AP s (t, \a\) = — / dv?e^AP(t,a) (139) 
2ir Jo 



If T is not large enough, then the peaks in the Fourier integral contain non-negligible contri- 
butions of the tails of the corresponding sine functions. In this case, the coefficients can be 
calculated from a set of linear equations obtained from eq. (134) for different times. 

55 For k= 1 very precise measurements of the Rabi oscillations have been performed recently (Brune, 
Schmidt-Kaler, Maali, Dreyer, Hagley, Raimond and Haroche [1996]), the peaked structure of the 
Fourier-transformed data being interpreted as a direct experimental verification of field quantization 
in a cavity. 

56 This condition is not satisfied, e.g., for even and odd coherent states as typical examples of 
Schrodinger-cat-\ike states. For even and odd coherent states, see Dodonov, Malkin and Man'ko, 
V.I., [1974]. For a review of Schrodinger cats, see Buzek and Knight [1995]. 
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(s = 0, 1, 2, . . .), it can be shown that they are related to the density-matrix elements 
by equations of the form of 



(for explicit expressions for Y£(t, \cx\), see Bodendorf, Antesberger, Kim, M.S., and 
Walther [1998]). Inverting eq. (140), which resembles, in a sense, eq. (112) in § 3.3.2., 
for each value of s then yields the density matrix of the cavity mode. Similarly to 
eq. (112), the inversion can be carried out numerically; e.g., by means of least-squares 
inversion (Appendix D.), choosing an appropriate set of values of |a|. 

In the two-mode nonlinear atomic homodyne detection scheme (Wilkens and Mey- 
stre [1991]) it is assumed that the signal cavity mode is mixed with a local-oscillator 
cavity mode according to the interaction Hamiltonian 



Equation (138) can then be found treating the effect of the local oscillator semiclas- 
sically (i.e., replacing the operator a L with a c number a^, a^ — >oll). In this case the 
scheme is obviously equivalent to an initial displacement a = —a^ of the density oper- 
ator of the cavity mode, so that the atomic-state inversion AP(t) is exactly given by 
eq. (138). In particular, when |a£,| is sufficiently large, then AP(t) can be rewritten 
as 



lated directly to the characteristic function of the Wigner function W{(f) of 

the cavity mode. Varying the interaction time and the phase ip L of a^, the whole 
function can be scanned in principle. Knowing $(/?), the Wigner function can 
then be obtained by Fourier transformation. 57 Since for appropriately chosen argu- 
ments $(/9) is nothing but the characteristic function of the quadrature-component 
distribution p(x, •&) for all values of d within a n interval (see Appendix B.5.), the 
density matrix in both a quadrature-component basis and the photon number basis 
can be reconstructed straightforwardly from eq. (142) (§§ 3.2. and 3.3.1.). Later it 
was found that the semiclassical treatment of the local oscillator restricts the time 
scale to times less than a vacuum Rabi period, because of the quantum fluctuations in 
the local-oscillator cavity mode (Zaugg, Wilkens and Meystre [1993]; Dutra, Knight 
and Moya-Cessa [1993]), and it was shown that this difficulty can be overcome when 
the atoms are weakly coupled to the local oscillator but strongly coupled to the signal 
(Dutra and Knight [1994]). 

A priori knowledge on the quantum state to be measured is required in mag- 
netic tomography (Walser, Cirac and Zoller [1996]). Here, the idea of quantum state 
mapping between multilevel atoms and cavity modes prepared in truncated states 
(Parkins, Marte, Zoller, Carnal and Kimble [1995]) is combined with a tomography 
of atomic angular momentum states by Stern-Gerlach measurements (Newton and 
Young [1968]). It is assumed that an angular-momentum degenerate two-level atom 



oo 




(140) 



n 



H' = tin a + (d + cll) + (o T + a L )a. 



(141) 




57 Recall that in the dispersive regime the Wigner function can be measured directly (Lutterbach 
and Davidovich [1997]; see § 2.4.), without any reconstruction algorithm. 
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passes adiabatically through the spatial profile of a classical laser beam [Rabi fre- 
quency: fl(t)} and, with a spatiotemporal displacement r > 0, through the profile 
of a quantized cavity mode [atom-cavity coupling: g(t — r)] such that the coupled 
atom-cavity system evolves according to the time-dependent Hamiltonian 

Je 

H(t) = hua^a + hcu eg | J e , m e ) ( J e , m e \ 

m e = ~J e 

-ihn(t) (e lULt A 1 - A[e- WLt ) + ihg(t-r) (a^A -Ala), (143) 

where the atomic deexcitation operators A , A±i are defined by 

A a = Yl \J 9 ,rn g ){J e ,me\Cl^X, (144) 

\m g \<J g , m e |< J e 

Ca,mg J m e being Clebsch-Gordan coefficients. 

It can be shown that if the time-dependent change of the Hamiltonian during the 
total interaction time is much less than the characteristic transition energies and if 
the delay and shape of the pulse sequences are chosen such that 

"^T* g(t-T)/n(t) oo, (145) 

then a coupled atom-cavity-field density operator qaf that can be factorized initially 
into a pure atomic state and a field state containing less than 2J g photons will be 
mapped to a product of atomic ground state superpositions and the cavity vacuum 

\J g , Jg-l)(Jg, Jg-1\® Q F ^ Q AF {t) QA ® |0)(0|, (146) 

with 

2J g -l 

Qf = (QF)mn\m)(n\ (147) 

m,n=0 

and 58 

2J g -l 

Qa= s m s n (g F )mn\Jg,Jg-m-l){Jg,J g -n-l\, (148) 

m,n=0 

s m , s n being possible sign changes. Hence, the original cavity-field density matrix 
(m\gF\n) is known if the final atomic density matrix (J g , J g — m — 1\qa\J 9 , J g — n — 1) 
is known. Atomic states of this type can be determined from magnetic dipole mea- 
surements using conventional Stern-Gerlach techniques (see § 4.5). 

3.6.2. Atomic beam deflection 

When a cavity mode is known to be in a pure state, then the expansion coefficients 
of the state in the photon-number basis can be inferred from the measured deflection 
of two-level probe atoms during their passage through the cavity (Freyberger and 
Herkommer [1994]; for an application of related schemes to the reconstruction of the 
transverse motional quantum state of two-level atoms, see § 4.4.1.). A narrow slit 



58 Note that with reverse adiabatic passage, an internal atomic state is preapred uniquely by 
reading out the cavity state. 
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put in front of one node of the standing wave (see Fig. 8) transmits the atoms only 
in a small region Ax <C A centered around x = (x, atomic center-of-mass position; 
A, wavelength of the mode). The dependence on x of the coupling strength k = k(x) 
— K'sm{kx) in eq. (65) can then be approximated by K — K'kx, k = 2ir/\, so that the 
interaction Hamiltonian of the combined system reads as 

H' = hn'kx (<r + a + a f (T_) . (149) 

It is further assumed that the system is initially prepared in a state = \^ f )\^a), 
where 

oo - 

\*f) = E c »' \*a) = 2- 1/2 / dx^(x) \x) (\g) + ce^e)) . (150) 

n=0 J 

The subscripts F and A label the states of the cavity field and the atom, respectively, 
ip(x) being the initial wave function of the atomic center-of-mass motion. Assuming 
c = 1 and calculating the temporal evolution of the state vector in the Raman-Nath 
regime, 59 it can be shown that the atomic momentum probability density for the 
interaction time t is given by 

oo 

P(p,t) = \hkJ2\±(P+V^hkK't)\ 2 \c n + e^c r ^\ 2 

n=l 

oo 

+lhkJ2\t(P-^ hkK ' t )\ 2 \ c n-e i ' p c^ 1 \ 2 + §M#(p)| 2 |c | 2 , (151) 

n=l 

where ip(p) is the initial wave function of the atomic center-of-mass motion in the 
momentum basis. For an initial position distribution |-?/>(x)| 2 , such as a Gaussian, 
with Ax <C A, and sufficiently strong interaction, i.e., n't — > oo, the corresponding 
shifted momentum distribution \ip(p± y^nhkn't)] 2 becomes sharply peaked at p = =F 
y/nhkn't, n= 1,2, .... In this case, (^(p)! 2 an d \ip(p ± y^nhkn't)] 2 obviously select 
in the measured momentum distribution P(p,t) a sequence of peaks whose heights 
are proportional to |c | 2 and \c n ± e l</3 c n __i| 2 , respectively, for n = and n> 0. In this 
way | c 1 2 and - if the measurement is performed at two phases ip - the sequence of 
products CqCi, cjc 2 , . . . can be obtained, from which (after truncating the state at a 
sufficiently large photon number n max ) the expansion coefficients c n can be calculated 
easily (up to an unimportant overall phase factor). 60 

To allow the cavity mode to be prepared in a mixed state, a running reference 
field aligned along the y-axis (i.e., perpendicular with respect to the standig-wave 
cavity field aligned along the x-axis and the z-direction of the moving atoms) can be 
included in the scheme such that the atoms initially prepared in the electronic ground 
state simultaneously interact with the two fields (Schneider, Herkommer, Leonhardt 
and Schleich [1997]). In the dispersive regime the effective interaction Hamiltonian 

59 In this approximation it is assumed that the transverse displacement of the atoms during the 
interaction time t is small compared to the wavelength of the field, so that the kinetic energy of 
the atomic center-of-mass motion can be disregarded in the Hamiltonian of the combined system. 
For times longer than the interaction time the now deflected atoms move freely, and the spatial 
distribution of the atoms on a screen put up far away from the resonator is a picture of their 
momentum distribution at the moment when they leave the interaction zone. 

60 For an analysis of a diffuse peak structure, see Freyberger and Herkommer [1994]. 
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for the system can then [on extending eq. (66) to two modes] be given by 

ii =n— — ( a_a + a*a — aa'a + a_ I + a — I (j-cr+a^ai, — a^a^(T + (T_ 1 

k(x)k l 



+h \(a_a + - a + a_) ( a + e^a^a)] . (152) 



5w 

Here the index L refers to the running reference field, which has the same frequency 
as the cavity field and can play a similar role as the local oscillator in a homodyne 
detector. In the calculations it is assumed that the atoms are initially prepared 
in a state \^a) of the form given in eq. (150) with c — 0, and the reference field is 
prepared in a strong coherent state and can be treated classically (o^— > «l = \aL\e llfL , 
\ul \ ~^ oo). 61 The cavity field can be prepared in an arbitrary mixed state described 
by a density operator g. Again making the approximation k(x) ~ n'kx and asuming 
that the initial atomic wave function ip(x,y) factorizes into a very narrow y-wave 
function = S(y — yo) and a Gaussian in the x-direction, it can be shown that the 
transverse momentum distribution P(p,t) of the atoms for an interaction time t is a 
(scaled) smeared quadrature-component distribution p(x, ip; rj) of the cavity mode, 

PM= ^kA x= ^kr™)' (153) 

R, = D,LK f hk\aL\t/5. Here, the quantum efficiency i] with which the quadrature- 
component distribution is measured is given by 

'=( 1 + 2Sy (154) 

and the phase <p is determined by the reference-field phase and the position y , ip — 
—ipL — ky . The result proves that the atomic probe carries the (smeared) quadrature- 
component information on the cavity field that corresponds to the reference angle 
ip. Hence, after having measured the transverse momentum probability distribution 
P(p,t) of the deflected atoms for a sequence of phases ^ of the (local-oscillator) 
reference field, the phase-space function P(q,p; s = 1 — i]^ 1 ) of the cavity mode can 
be reconstructed from eq. (75) by means of tomographic imaging (§ 3.1.). Note that 
the larger the effective coupling constant k is, the closer is P(p, q; s) to the Wigner 
function. However the initial position uncertainty must not exceed a fraction of the 
wavelength of the cavity mode, because of the approximations made. Needless to say, 
the interaction time t must be small compared to the decay time of the cavity mode. 



3.7. Alternative proposals 

Since the quantum state of a radiation field mode is known when its density-matrix 
elements g mn in the photon-number basis are known, it was proposed to measure them 
directly (Steuernagel and Vaccaro [1995]). The diagonal elements can be measured 
by direct photodetection in principle. In order to measure the off-diagonal elements, 

61 Note that an atom that is initially in the electronic ground state (c = 0) stays in the electronic 
ground state, because of the effective Hamiltonian (152). 
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it was proposed to prepare a probe field in a superposition of two photon-number 

states, 

\orrm) = N a (\m) + a\n)), (155) 

N a = (1 + |a| 2 )~ 1//2 , and combine it with the signal mode at a beam splitter. The 
observed joint photon-number probabiltity in the two output channels of the beam 
splitter then reads as (m > n) 

Pq,k+m-q ~ (&fcj|£|&fcj), 

/ — k — m — n, where \bki) is again a state of the type given in eq. (155), so that 

(hi\Q\hi) = Nl [^o kk + \b\ 2 Q U + bQ H + b*Qu} . (157) 

When for chosen difference k — I = n — m the diagonal elements Qu and ^ are 
known from a direct photon- number measurement, then the off-diagonal elements 
Qkk+m-n, k— 0,1,2,..., can be obtained from two (ensemble) measurements for 
different superposition parameters. The whole density matrix can then be obtained 
by means of a succession of (ensemble) measurements, varying the difference m — n 
of photon numbers in the reference state (155) from measurement to measurement. 
The presently unresolved problem, however, consists in building an apparatus that 
prepares travelling waves in superpositions of two photon- number states \m) and \n) 
for arbitrary difference m — n in a, controllable way. 

For a radiation- field mode that is known to be prepared in a pure state \^/) and 
that contains only a finite number of photons, eq. (107), it was proposed to extend 
the method of photon chopping outlined in § 2.1.7. to a complete determination of 
the expansion coefficients c n (Paul, H., Torma, Kiss and Jex [1996a, 1997a]). 62 First 
the photon-number distribution p n = \c n \ 2 is determined according to eqs. (52) - (54). 
Then the balanced 2iV-port apparatus (N > n max ) is used for mixing the signal mode 
with N — 2 reference modes in the vacuum state and one reference mode prepared in 
a coherent state \a) (in place of the N vacuum reference inputs in § 2.1.7.). Using the 
selected statistics when only photons in the first N/2 output channels are observed, 
it can be shown that the probability p n (N) of detecting n photons with at most one 
photon in each channel is given by (Paul, H., Torma, Kiss and Jex [1996a]) 

2 

(158) 

Determining the coincidence probability distribution p n (N) for two coherent states 
with different phases, the phases ip n of the expansion coefficients c n = \c n \e t(Pn can be 
calculated step by step (for known amplitudes \c n \ and up to an overall unimportant 
phase) from the two measured coincidence-event distributions. More involved formu- 
las are obtained in the case when the complete coincidence statistics is included in the 
analysis; i.e., photons are allowed to be observed in any output channel and realistic 
photodetection is considered (Paul, H., Torma, Kiss and Jex [1997a]). It turns out 
that the reconstruction scheme with two coherent states requires photodetectors that 

62 For a proposal to extend the method to a two- mode signal field, see Paul, H., Torma, Kiss and 
Jex [1997b]. 



(156) 



Pn(N) = 



'N/2\ e-w 
n I N n 



n 



E r)y/(n-m)\orc n . 

m=0 \ m / 
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can discriminate between 0, 1, 2, . . . photons, i.e., detectors that have rather low quan- 
tum efficiency. Using avalanche photodiodes, which can have high quantum efficiency, 
the reconstruction might become rather involved, because additional reference beams 
with different phase properties must be used. Similarly, when a mixed quantum state 
is tried to be reconstructed, then extra reference beams (reference phases) must be 
used, the maximum number of reference beams being limited by the cutoff in the 
photon number. 

It can also be shown (Bialynicka-Birula and Bialynicki-Birula [1994]; Vaccaro and 
Barnett [1995]) that when a radiation-field mode is prepared in a pure state |\&) that is 
a finite superposition of photon-number states, eq. (107), then it can be reconstructed 
from the photon-number distribution p n = \c n \ 2 and the (Pegg-Barnett) 63 truncated 
canonical phase distribution Ppb{4>) — \^{4>) | 2 > with 64 

1 ^max 

= \m\e tm = —j E \cn\^ e**- (159) 

V ""max ~r J- n= Q 

Taking eq. (159) at 2(n max + 1) values of (f>, such as <pi = In/ ( n max + 1), the resulting 
equations can be regarded as conditional equations for the unknown phases ip n , pro- 
vided that all absolute values and \c n \ are known. Apart from the fact that 
the quantum state must be known to be a pure state, the question remains of how to 
obtain the phase statistics (see also § 3.8.3.). 

3.8. Reconstruction of specific quantities 

Since the density matrix in any basis contains the full information about the quantum 
state of the system under consideration, all quantum-statistical properties can be 
inferred from it. Let F be an operator whose expectation value, 

Fnm Qmn i 

(160) 

mn 

is desired to be determined. One may be tempted to calculate it from the recon- 
structed density matrix (or another measurable quantum-state representation). How- 
ever, an experimentally determined density matrix always suffers from various inac- 
curacies which can propagate (and increase) in the calculation process (cf. § 3.9.1.). 
Therefore it may be advantageous to determine directly the quantities of interest from 
the measured data, without reconstructing the whole quantum state. In particular, 
substituting in eq. (160) for g mn the integral representation (91) one can try to obtain 
an integral representation 

(F) = dip J dxK F (x,<p)p(x,<p), (161) 

with 

K F (x, v?) = £ F nm f mn (x) e*( m -^, (162) 

mn 

suited for direct sampling of (F) from the quadrature component distributions p(x, ip) 
[provided that both the sampling function K F (x, ip) and the integral in eq. (161) exist]. 

63 See Pegg and Barnett [1988]. 

64 Note that ip(4>) = ($pb(</>)| , I'), where |^} is given by eq. (107), and the truncated phase state 
|*pb(0)) reads as |$ PS (0)) = (rw + 1)" 1 / 2 £^= a x e - in *|n) (cf. § 3.8.3.). 
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3.8.1. Normally ordered photonic moments 

It is often sufficient to know some moments of the photon creation and destruction 
operators of a radiation-field mode rather than its overall quantum state. Let us 
again consider the measurement of the quadrature-component probability distribu- 
tions p(x, if) in balanced homodyning (§ 2.1.2.). Then the normally ordered moments 
of the creation and destruction operators, {o) n a m ) can be related to p(x, </?) as 65 

(a*dT) = 7r2(w+mV2(w+m) , I ^ J d*p(s, <f) M^, (163) 

where H n (x) is the Hermite polynomial (Richter, Th., [1996b]). 66 Equation (163) 
is just of the type given in eq. (161) and offers the possibility of direct sampling of 
normally ordered moments (a^ n a m ) from the homodyne data. It is worth noting that 
knowledge of p(x, <p) at all phases within a n interval is not necessary to reconstruct 
{a) n a m ) exactly, and therefore the (p integral in eq. (163) can be replaced with a sum. 
It was shown that any normally ordered moment (a^ n a m ) can already be obtained 
from p(x, ip) at N = n + m + l discrete different phases ipk (Wiinsche [1996b, 1997]). 

The method is especially useful, e.g., for a determination of the moments of photon 
number. Note that for finding the mean number of photons (n) from the Fock-basis 
density matrix a relatively large number of measured diagonal elements g nn must be 
included into the calculation, in general, each of which being determined with some 
error. After calculation of the sum (n) = J2 n n Qnn the error of the result can be too 
large to be acceptable, and for higher-order moments severe error amplification may 
be expected. Equation (163) can easily be extended to nonperfect detection (rj < 1), 
since replacing p(x, ip) with p(x, ip; 77) simply yields ^( n+m )/ 2 ^ n d m ) , 67 

Provided that the whole manifold of moments (a) n a m ) has been determined, then 
the quantum state is known in principle (Wiinsche [1990, 1996b]; Lee [1992]; Herzog 
[1996b]). 68 To be more specific, the density operator can be expanded as [cf. eq. (69)] 

00 

q = ]T a k>l (a^ k a l ), (164) 

k,l=0 



where 



(-W-j){k-j\ ( . 

a k,i = 2^ 1 ( 165 ) 



3 



^ j\y/(k-j)\(l-j)\ 



[{k, 1} = min (k, /)], provided that (a) k a l ) exists for all values of k and / and the series 
(164) (in chosen basis) exists as well. Since all the moments do not necessarily exist for 



65 Equation (163) can be proven correct if both sides are expressed in terms of the density matrix 
in the photon-number basis and standard summation rules for the Hermite polynomials are used. 

66 A more involved transformation was suggested by Bialynicka-Birula and Bialynicki-Birula [1995]; 
see footnote 77. 

67 Recall that an imperfect detector can be regarded as a perfect detector with a beam splitter in 
front of it, so that the destruction operator a of the mode that is originally desired to be detected 
is transformed according to a beam-splitter transformation (§ 2.1.1.), a{rf) = y/rja + y/1 — r\ b, where 

b is the photon destruction operator of a reference mode prepared in the vacuum state. 

68 This is also true for other than normally ordered moments, such as symmetrically ordered 
moments (Band and Park [1979]; Park, Band and Yourgrau [1980]). 
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any quantum state, and the series need not necessarily converge for existing moments, 
the quantum-state description in terms of density matrices is more universal than that 
in terms of normally ordered moments. 69 

The extension of the method to the reconstruction of normally-ordered moments 
of multimode fields from joint quadrature-component distributions is straightforward. 
Moreover, normally-ordered moments of multimode fields can also be reconstructed 
from combined distributions considered in § 2.1.3. (Opatrny, Welsch and Vogel, W., 
[1996, 1997a]; McAlister and Raymer [1997a,b]; Richter, Th., [1997a]). For simplicity, 
let us restrict attention here to two-mode moments (a| rai a} 2 2 a™ 1 a™ 2 ) (the subscripts 
1 and 2 refer to the two modes) and assume that the weighted sum x — Xi(ipi) cosa-l- 
x 2 (^2) si n ex, eq. (29), is measured in a homodyne detection scheme outlined in § 2.1.3. 
(see Fig. 4). Measuring the moment (x n ) for n+ 1 values of a, a set of n + 1 linear 
algebraic equations can be obtained whose solution yields the two-mode quadrature- 
component moments (xi~ k {(pi)x\{(p 2 )) (k — 0, 1, 2, . . . , n). Varying (fi and </? 2 , the 
procedure can be repeated to obtain the dependence on tpi and ip 2 of the two-mode 
quadrature-component moments. These can then be used to determine the normally 
ordered moments of the photon creation and destruction operators of the two modes 
step by step. In a closed form, (a{ ni a 2 n2 a™ 1 a™ 2 ) can be obtained from the combined 
quadrature-component distribution ps(x, a, ipi, <p 2 ) as 



h 



p s {x : a : (p 1 ,cp 2 ) 



/r f"K f"K 

dx da dy?i / dip 2 
Jn Jo Jo 

x H ni ^^ 2 {x) F^^^ia) e ^-^) e i(n 2 -m^ {m) 

[C nm = n\m\ / (ir2( n+m ^ 2 (n + m)!)]. The functions F^\a) form a biorthonormal system 
to the functions G^\at) = (f) cos'~ fe a sin fc o; in some interval, so that J n da F^\a>) 



r fc ~ \ k J 

G$(a) = 5kk'- Note that since for given / there is a finite number of functions G^)(at), 
the system of functions F^\a) is not uniquely determined and can be chosen in 
different ways. In particular when only the phase difference Aip = ip 2 — y?i is controlled 
and the overall phase </?i is averaged out (this is the case when the signal and the local 
oscillator stem from different sources), then the moments (a\ ni a^' 1 'a™ 1 a™ 2 ) can be 
reconstructed for n\ — mi = m 2 — n 2 . If both phases ipi and ip 2 are averaged out, then 
those moments (a| ni a} 2 n2 ^a" 1 a 2 2 ) can still be obtained, which carry the information 
about the photon-number correlation in the two modes. 

The method was used to experimentally demonstrate the determination of the 
ultrafast two-time photon number correlation of a 4 ns optical pulse (McAlister and 
Raymer [1997a]). In the experiment, the local-oscillator pulses are derived from 
a Ti:saphire-based laser system which generates ultrashort, near transform-limited 
pulses (150 fs) at a wavelength of 830 nm and a repetition rate of 1 kHz. The signal 
is from a single-spatial-mode superluminescent diode. The broadband emission at 
830 nm is spectrally filtered to produce a 4 ns pulse having a 0.22 nm spectral width. 



69 An example of noncxisting normally ordered moments is realized by a quantum state whose 
photon-number distribution behaves like p n ^n~ 3 for n^oo. Even though it is a normalizable state 
with finite energy, its moments (h k ) do not exist for k > 2. For a thermal state the relation (a^ n a m ) 
= 5 nm n\n n is valid (n, mean photon number), which implies that the series (164) for g mn does not 
exist. Note that the problem of nonconvergence of the series (164) may be overcome by analytic 
continuation. 
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Figure 15: The second-order coherence, eq. (167), experimentally determined 
via balanced four-port homodyne detection (dots) and from the measured 
optical spectrum (solid line). The value of t\ is set to occur near the maximum 
of the signal pulse. (After McAlister and Raymer [1997a].) 



For each value of the relative delay At = t\ — 1 2 between two local-oscillator pulses the 
phase-averaged combined quadrature-component distribution ps(x, a) = (2n)~ 2 J dipt 
J dip2Ps(x, a, <fi, (£2) is measured for three different values of a: 0, 7r/2 and 7r/4. From 
these the normalized second-order coherence function 



is computed (Si and a 2 being the photon destruction operators of the spatial-temporal 
modes defined by the local-oscillator pulses centered at times t\ and t 2 \ This 
experiment represents an extension of measurements of the Hanbury-Brown-Twiss 
correlations 70 to a sub-picosecond region. Results are shown in Fig. 15. 



3.8.2. Quantities admitting normal-order expansion 

The basic relation (163) can also be used [similar to eqs. (160) - (162)] to find sampling 
formulas for the mean values of quantities that can be given in terms of normally 
ordered moments of photon creation and destruction operators (D'Ariano [1997b]). 
Let us consider an operator F and assume that it can be given by a normal-order 
expansion as 



9 




(a\a 2 aia2) 



(167) 



(a\ai) (a\a 2 ) ' 



CO 



F — ~y ' fnm 0^ & 



(168) 



n,m=0 



70 See Hanbury Brown and Twiss [1956a, b, 1957a, b]; for a more discussion of the Hanbury Brown 
and Twiss experiments, see Pefina [1985], Mandel and Wolf [1995]. 
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The mean value (F) can then be obtained from p(x, ip; rj) according to 71 

(F) = J dy? J dxK F (x,<p;rj)p(x,<p;rj), (169) 
where the integral kernel 



K F (x,ip;rj) = Tr FK(x,<p;ri) 



(170) 



[K(x, ip; rj) being defined in eq. (97)] can be given by 

K F (x,<p;r,) = —7^=f J dye-^/ 2 G F (y+z^/2^x,^ , (171) 



with 



G F (z,<p) 




/ j J nm 
n,m=0 



-iz) 



n+m ^i(m—n)tp 



(172) 



Here it is assumed that the integral kernel K F (x,ip;r]) exists and the x integral in 
eq. (169) exists as well. 72 . The integral kernel (171) exists if for z^oo the function 
G F (z,ip) grows slower than e vz / 2 , so that the integral in eq. (171) converges (for 
examples, see D'Ariano [1997b]; D'Ariano and Paris [1997b]). 73 The integral kernel 
(171) is therefore applicable to such quantum states whose (smeared) quadrature- 
component distributions p(x, ip; rj) tend to zero sufficiently fast as \x\ goes to infinity 
such that eq. (169) converges even if K F (x,ip;r]) increases with \x\. At this point it 
should be noted that K F (x,ip;rj) is determined only up to a function Q(x,ip); i.e., 
K F (x, ip; rj) in eq. (169) can by replaced with K' F (x, ip; rj) = K F (x, ip; rj) + Q(x, ip) such 
that w 

J dip J dx Q(x, ip) p(x, ip; rj) = (173) 

for any normalizable quantum state. Hence, if the integral kernel K F (x,ip;rj) that is 
obtained from eq. (171) is unbounded for \x\ ^oo such that the x integral in eq. (169) 
does not exist for any normalizable quantum state, it cannot be concluded that F 
cannot be sampled from the quadrature-component distributions of any normalizable 
quantum state, since a different, bounded kernel may exist. 



3.8.3. Canonical phase statistics 

The quantum-mechanical description of the phase and its measurement has turned 
out to be troublesome and is still a matter of discussion. Many papers have dealt 
with the problem and an extensive literature is available (for reviews, see Luks and 

71 Recall that when the c- number function F(a;s) that is associated with an operator F in s order 
exists for s < — 1, then (F) can be obtained from the phase-space function P(a; s) (measurable, e.g., 
in balanced eight-port homodyning; see § 2.1.4.) according to eq. (B.22), provided that the integral 
exists. 

72 Equations (169) - (172) can be obtained by taking the average of eq. (168), substituting for 
(a'"a m ) the right-hand side (multiplied by ry _ (™+ m )/ 2 ) of eq. (163), and using the generating function 
of the Hermite polynomials with the argument (i/y/^rj)(d/dz). 

73 Since the integral kernel Kp(x,ip;rj) does not depend on ip if F is a function of the photon- 
number operator, this case can be regarded as the first realization of Hclstrom's quantum roulette 
wheel (D'Ariano and Paris [1997b]). 
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Perinova [1994], Lynch [1995], Royer [1996], Pegg and Barnett [1997]). Here we con- 
fine ourselves to the canonical phase that is obtained in the attempt - in close analogy 
to the classical description - to decompose the photon destruction operator into am- 
plitude and phase such that a = E\fh. The (nonorthogonal and unnormalizable) 
phase states, 

|0> = fV"V> (174) 

n=0 

(London [1926, 1927]), which are the right-hand eigenstates of the one-sided unitary 
operator E, E\4>) = e l<t> \4>) , can then be used to define - in the sense of a POM - 
the canonical phase distribution P(0) = (11(0)), with 11(0) = |0)(0|. 74 Note that the 
Hermitian operators C— (E+E J( )/2 and S = (E — &)/(2i) are the Susskind-Glogower 
sine and cosine operators (Susskind and Glogower [1964]). 75 

Various proposals have been made to measure P(4>) using homodyne detection or 
related schemes. 76 So the Wigner function or the density matrix reconstructed from 
the homodyne data could be used to infer the phase statistics in principle (Beck, 
Smithey and Raymer [1993]; Smithey, Beck, Cooper and Raymer [1993]; Breitenbach 
and Schiller [1997]). Further it was proposed to obtain P(0) = (11(0)) from a homo- 
dyne measurement of (all) normally ordered moments (aJ m+n a n ) , because 11(0) has a 
series expansion of the type given in eq. (168), 



n(0) = (2tt) 



-1 



oo oo 

1 + H (e im *a tm+n a n + e- ir ^a tn a m+; ' 

m=l n=0 



(175) 



(Bialynicka-Birula and Bialynicki-Birula [1995]). 77 Unfortunately (11(0)) cannot be 
given (even for r] — 1) by an integral relation of the form (169) suited to statistical 
sampling, because the integral kernel does not exist (see footnote 79). It was therefore 
proposed to introduce parametrized phase distributions P(0, e) such that for any e > 
P(0, e) can be obtained by direct sampling from p(x, ip), and P(0, e) — >P(0) for e^O 
(Dakna, Knoll and Welsch [1997a,b]). The sampling functions can then be obtained 
as convergent sums of the sampling functions for the density-matrix elements in the 
Fock-basis (see § 3.3.). Since the value of the parameter e for which P(0, e) f»P(0) is 
determined by the number of photons at which the quantum state under study can 
effectively be truncated, the method is state-dependent. Further, it was proposed to 
measure the canonical phase by projection synthesis, mixing the signal mode with 
a reference mode that is prepared in a quantum state such that, for appropriately 
chosen parameters, the joint-photon-number probabilities in the two output channels 
of the beam splitter directly correspond to the canonical phase statistics of the signal 
mode (Barnett and Pegg [1996]; Pegg, Barnett and Phillips [1997]). Apart from the 



74 For a two-mode orthogonal projector realization in the relative-photon-number basis, see Ban 
[1991a-d, 1992, 1993]; Hradil [1993]. 

75 For constructing the operators corresponding to classical phase-dependent quantities, see Bergou 
and Englert [1991]. 

76 The proposal was made to regard the quadrature-component distribution at x = as (unnormal- 
ized) phase distribution, P{4>) ~p(a; = 0, 4>) (Vogcl, W., and Schleich [1991]), which yields a phase 
statistics that is quite different from the canonical phase statistics in general (for an improvement, 
see Buzek and Hillery [1996]). 

77 Here it was proposed to obtain (a) m+n a n ) from the homodyne data via the (2n + m)th derivative 
of J Q 27r d^e- m ^e- A2 / 2 (exp[iAx(^-7r/2)]). 
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direct photon-number measurement needed and the fact that the quantum state under 
study must again be truncated, so that the method is state-dependent, the difficult 
problem remains to design an apparatus that produces the reciprocal binomial states 
needed. 

The problem of homodyne measurement of the canonical phase can be solved 
when the exponential phase moments fy k [i-e., the Fourier components of P{4>)\ are 
considered and not the phase distribution itself, 



+00 



P(0) = (27T)- 1 £ e~ ik H k , 

k=—oo 



(176) 



where ^ k = (E k ) if k > 0, and ^ k = V*_ k if k < (tf = !)■ Then it can be shown that 
(k>0) 

^k = j Q d<p J dxe ik *K k (x)p(x,<p) (177) 

(Opatrny, Dakna and Welsch [1997, 1998]; Dakna, Opatrny and Welsch [1998]). The 
integral kernel K k (x) (Fig. 16) can be used for sampling the exponential phase mo- 
ments from the homodyne data for any normalizable state. 78 In particular, K k (x) 
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Figure 16: The x-dependent part K k (x) of the sampling function K k (x,ip) 
= e lkip K k {x) for the determination of the exponential moments ^ k = (E k ) 
of the canonical phase from the quadrature-component distributions p(x, tp) 
according to eq. (177) for various (a) odd and (b) even k. (After Dakna, 
Opatrny and Welsch [1998].) 



rapidly approaches the classical limit 79 

r,n f i(-l) (fe - 1)/2 A;signa; if k odd, 

[ (27r) L (-iy k+2 >' 2 k\nx if k even 

as I a; I increases, and it differs from the classical limit only in a small interval around 
the origin. It is worth noting that the method applies to quantum and classical 



78 Note that when K k {x) is calculated according to eq. (171), then it docs not apply to all normal- 
izable states (for analytical expressions and properties of Kk (x) , including nonperfect detection, sec 
Dakna, Opatrny and Welsch [1998]). 

79 Already from the classical kernel (178) it is seen that J2t=-oo e lk ^'^K k {x) does not exist, and 
hence P{4>) cannot be obtained from p(x, tp) by means of an integral transformation of the form 
given in eq. (177). 
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systems in a unified way and bridges the gap between quantum and classical phase. In 
particular, the integral kernel K^ l \x) as given in eq. (178) is nothing but the integral 
kernel for determining the radially integrated Wigner function, which reveals that in 
the classical limit the canonical phase distribution is simply the radially integrated 
probability density for the complex amplitude a. Note that any radially integrated 
propensity prob(o;), such as the radially integrated Q function measurable, e.g., in 
balanced eight-port homodyne detection, can be used to define an operational phase 
probability density (Noh, Fougeres and Mandel [1991, 1992a,b, 1993a,b,c]), which 
in the quantum regime differs from the canonical phase distribution in general (for 
further readings, see Turski [1972]; Paul, H., [1974]; Shapiro, J.H., and Wagner [1984]; 
Hradil [1992, 1993]; Vogel, W., and Welsch [1994]; Leonhardt, Vaccaro, Bohmer and 
Paul, H., [1995]; Leonhardt [1997c]). 

3.8.4. Hamiltonian and Liouvillian 

So far, measurement and reconstruction of quantum-state representations and aver- 
ages of particular quantities at certain time have been considered. The quantum state 
of an object at chosen time t is of course a result of state evolution from an initial time 
to, the state evolution being governed by the Hamiltonian H of the system or a Liou- 
villian L in the more general case of the system being open. Since the Liouvillian of 
an object expresses in most concentrated form the dynamics of the object, knowledge 
of the Liouvillian is essential for understanding the behaviour of the object, and the 
question may arise of how to experimentally determine its form. To answer the ques- 
tion, it was proposed to appropriately apply quantum-state reconstruction routines, 
such as direct sampling of the density matrix in balanced homodyning (D'Ariano and 
Maccone [1997]). 

Let us consider a system that is initially prepared in some known state g- in = g(to) 
and assume that in the further course of time it evolves to a state g ont = g(t) at time t, 

Pont = G g- m , (179) 

where for a system that is homogeneous in time the superoperator G = G(t, to) can 
be given by the exponential of the (time-independent) Liouvillian L of the system, 

G = e L{t - ta) . (180) 

The superoperators G and L can then be obtained by measuring (reconstructing) 
g ont for various probe inputs f3; n . A typical example may be a radiation-field mode 
which (at time to) is fed into a nonlinear resonator-like equipment giving rise to 
amplification and damping. The quantum state of the outgoing mode (at time t) can 
then be related to that of the incoming mode according to eq. (179). In the Fock 
basis, eq. (179) reads as 

oo 

(m\g out \n) = £ G l * n {l\g in \k). (181) 

l,k=0 

Now let us assume that in a succession of measurements the quantum states of the 
outgoing mode are reconstructed for an appropriately chosen (overcomplete) set of 
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input quantum states g^, /i — 1,2, . . ., such that for chosen m and n the block of 
equations 

oo 

HeS» = E G£»<i|eU*> (182) 

Z,fc=0 

can be used as conditional equations for the matrix elements G l £ n . Repeating the 
procedure for all values of m and n, the G matrix and [by using eq. (180)] the L 
matrix can be calculated in principle. Apart from a truncation of the Hilbert space 
in numerical implementations, the main problem that remains in practice is of course 
the (controlled) initial-state preparation such that all (relevant) matrix elements of 
the Liouvillian can be probed. 

In particular, when the Liouvillian is phase-insensitive, then mixtures of Fock 
states evolves into mixtures of Fock states. For this case a setup was proposed which 
uses correlated twin beams produced by a nondegenerate parametric amplifier in 
combination with conditional photon-number measurement on one beam in order to 
prepare the other beam in (random) photon-number states. These states can then be 
used as inital states for probing the Liouvillian. To realize the scheme with available 
techniques, it should be noted that the reconstruction algorithm can be extented such 
that the quantum efficiency of the photon-number measurement need not be unity 
(D'Ariano [1997c]). 

3.9. Processing of smeared and incomplete data 

In practice, there are always experimental inaccuracies that limit the precision with 
which the quantum state in a chosen representation can be determined. Typical 
examples of inaccuracies 80 are data smearing owing to nonperfect detection, finite 
number of measurement events, and discretization of continuous parameters, such as 
the quadrature-component phase in balanced homodyning. In particular, the latter is 
an example of leaving out observables A, in the expansion (69) of the density operator 
q. Whereas in balanced homodyning the distance between neighboring phases can 
be diminished such that the systematic error is reduced, in principle, below any 
desired level, there are also cases in which a principally incomplete set of observables 
is available. Then either additional knowledge of the quantum state is necessary to 
compensate for the lack of observables or other principles must be used to reconstruct 
the density operator according to the actual observation level. 

3.9.1. Experimental inaccuracies 

Since in a realistic experiment the quantum efficiency r\ is always less than unity, the 
probability distributions of the measured quantities are always more or less smeared, 
because of losses. The problem of compensating for losses in direct sampling of the 
quantum state from the quadrature- component distributions measurable in balanced 
homodyning (§ 2.1.2.) has been studied widely. As mentioned in § 2.3., active loss 
compensation may be realized, in principle, by means of a squeezer, such as a de- 
generate parametric amplifier (Leonhardt and Paul, H., [1994a]; see also Leonhardt 

80 For limits upon the measurement of state vectors expressed in terms of channel capacities for the 
transmission of information by finite numbers of identical copies of state vectors, see Jones, K.R.W., 
[1994]. 
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and Paul, H., [1995]). When the signal mode is preamplified before detection by 
means of a squeezer [eq. (63) with e l ^d s in place of d\], then it can be shown that for 
appropriately phase matching the (scaled) quadrature-component distribution of the 
preamplified signal measured with quantum efficiency rj reads as 81 



(g, amplification factor), with p(x; rf) being given by eq. (20). From eq. (183) together 
with eq. (20) it is seen that r\ is effectively replaced with gr\j ( 1 — rj+gn), and hence ^fg 
p'iyfgx, if] rf) tends to p(x, ip) as the amplification becomes sufficiently strong. 82 Note 
that the degree of improving the quantum efficiency is limited by the realizable mode- 
matching in degenerate parametric amplification, which deteriorates with increasing 
pump strength. 

Without active manipulation of the signal, it has been well established that above 
a lower bound for rj, it is always possible to compensate for detection losses, introduc- 
ing a modified sampling function that depends on 77 such that performing the sampling 
algorithm on the really measurable (i.e., the smeared) quadrature-component distri- 
butions p(x, p; rf) yields the correct quantum state (D'Ariano, Leonhardt and Paul, 
H., [1995]). This bound depends on the chosen state representation in general. In 
particular, reconstruction of the density matrix in the photon-number basis is always 
possible if 77 > 1/2. In this case, 



f o d( pf dxK mn (x,p;r])p(x,p;r]), (184) 



where the now 77-dependent sampling function K mn (x, p; rf) is defined by eq. (96) 
[together with eq. (97)]. 

Another approach to the problem of loss compensation is that the sampling func- 
tion is left unchanged [i.e., K mn (x, p; rf) ->• K mn (x, p) = K mn (x, p;r) = l) in eq. (184)] 
and, at the first stage, the density matrix g m n(v) °f the quantum state that corre- 
sponds to the smeared quadrature-component distributions p(x, p; rf) is reconstructed. 
After that, at the second stage, the true density matrix g mn = g(r] = 1) is calculated 
from Qmniv) (Kiss, Herzog and Leonhardt [1995]; Herzog [1996a]). Extending the 
Bernoulli transformation (C.3) and the inverse Bernoulli transformation (C.4) to off- 



81 Note that squeezing the local oscillator has the same effect as antisqueezing the signal field under 
the assumption that the coherent component is large (Kim, M.S., and Sanders [1996]; see also Kim, 
M.S., [1997b]). 

82 For details, see also Vogel, W., and Welsch [1994], and for a discussion of eqs. (183) and (20) 
in terms of moments of the measured quadrature components and those of x(ip), see Marchiolli, 
Mizrahi and Dodonov [1997]). 
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diagonal density-matrix elements, g mn can be calculated from g m n(v) as83 

- = """' + " ,/2 £ C?f (if 2 ( l - f < 185 > 

It is worth noting that eq. (185) is exact. In other words, when (for precisely given 
overall detection efficiency) Q mn {rj) is exactly known, then g mn can always be calcu- 
lated precisely, irrespective of the value of rj. 

In practice however, r\ is not known precisely in general, and the experimentally 
determined values of g m n(v) always differ from the exact ones. In particular, when 
7] < 1/2 and the error of g mn (v) does not vanish with increasing m and n, and when 
there is no a priori information about the quantum state, such as the photon number 
at which it can be truncated, then g mn cannot be obtained - for chosen number 
of measurements - from the measured Qm(v), because of error explosion. A typical 
example of such an error is the statistical error with which g mn (r)) is sampled from the 
quadrature-component distributions measured in balanced homodyning. In this case 
loss compensation (for arbitrary quantum states) is only possible if rj > 1/2 (D'Ariano 
and Macchiavello [1998]). 84 Clearly when the quantum state to be reconstructed can 
be truncated at some maximum photon number n max , then the sum in eq. (185) is 
finite and the density- matrix elements g mn can also be obtained for rj < 1/2, the 
accuracy being determined by that of Q mn {rj). 

Limitations on variables in real experiments always give rise to systematic errors. 
This is the case, e.g., when the quadrature components x(ip) are measured at discrete 
phases or when only some part of the ir interval can be scanned experimentally 
(for an example, see § 4.4.1.). In particular, in balanced homodyning the quadrature- 
component distribution p(x, (p) is measured at a finite number of phases ipk within a ir 
interval and finite ^-resolution. When p(x, (p) is measured precisely at N equidistant 
phases ipk = (ir/N)k, then the reconstructed density- matrix elements are given by 

Q mn {N) = — ^ y dxK mn (x,ip k )p(x,ip k ), (186) 

in place of eq. (91). A measure of the systematic error is the difference Ag mn = 



83 Replacing in eq. (185) r\~ x with r\ yields g mn (j]) i n terms of g m n, which corresponds to a gen- 
eralization of the Bernoulli transformation (C.3) to off-diagonal density-matrix elements. Since 
an imperfect detector can be regarded as a perfect detector with a beam splitter in front of it, 
g(r]) and g can be related to each other applying the beam-splitter transformation (7), g(rj) = 
2 (0|e^ 2 f3e-^ £2 |0) 2 , with L 2 according to eq. (6) and cos 2 (/3/2) = 77 (the mode indices 1 and 2, 
respectively, refer to the signal and vacuum inputs). It can then be proved easily that g(rj) = rj L g, 
where L= ^(a^ag + ga^a — 2aga { ), which corresponds to the time evolution of a damped harmonic 
oscillator at zero temperature (77— >e~*). 

84 When the density matrix at initial time t = of a signal mode that undergoes phase-insensitive 
damping or amplification is tried to be reconstructed from the quadrature-component distributions 
p(x, <p, t; 77) measured at time t > 0, then the additional (phase-insensitive) noise gives rise to a 
modified overall quantum efficiency 77*, so that 77* > 1/2 must be valid in order to compensate for 
the losses (D'Ariano [1997b]; D'Ariano and Sterpi [1997]). 
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Qmn(N) - g mn , which reads as 85 

oo 

Ag mn = J2J2 G T0ki, with k-l = m-n±2jN (187) 

3=1 k,l 

[G™ n = 7T / dx f m n(x) g k i{x)\ (Leonhardt and Munroe [1996]; Leonhardt [1997b]). 
Equation (187) reveals that when the quantum state can be truncated such that 

g mn = for \m-n\ > N, (188) 

then the density matrix-elements g mn for \m — n\ < N can be reconstructed precisely 
from p(x, ip) at N phases ip k according to eq. (186). In particular, one phase is required 
to reconstruct a completely dephased quantum state that contains only diagonal 
density-matrix elements. Since any quantum state can be approximated to any desired 
degree of accuracy by setting g mn m for m(n) > n max , there is always aniVft! n max 
for which the condition (188) can be assumed to be satisfied, so that all (essentially 
nonvanishing) density-matrix elements can be reconstructed precisely. If iV is not 
known a priori, the quantum state can be reconstructed in an iterative way, increasing 
N in a sequence of (ensemble) measurements until |A^ mn | is sufficiently small. Note 
that any normally ordered moment (a^ n a m ) can exactly be reconstructed from p(x, tp) 
at N — n + m + 1 phases (Wiinsche [1996b]; see § 3.8.1.). 

The methods for quantum-state reconstruction are based on ensemble measure- 
ments; i.e., on a sequence of individual measurements carried out on identically pre- 
pared systems. 86 Since the number of individual measurements is principally finite, 
the measured quantities are always estimates of the true ones. Hence all the quan- 
tities that can be inferred from the measured quantities are also estimates, and the 
statistical error with which the original quantities are measured propagates to the 
quantities inferred from them. In particular, when the quantities which are desired 
to be determined can be sampled directly from the measured data, then the statis- 
tical error can be estimated straightforwardly by also using the sampling method. 
The problem of statistical error in quantum-state reconstruction has been studied in 
a number of papers, with special emphasis on balanced four-port homodyne detec- 
tion (Leonhardt, Munroe, Kiss, Raymer and Richter, Th., [1996]; D'Ariano [1997b]; 
D'Ariano and Paris [1997a]; D'Ariano, Macchiavello and Sterpi [1997]; Leonhardt 
[1997c]). Let us consider a quantity F and assume that it can be sampled directly 
from the quadrature-component statistics according to eq. (161). When in an exper- 
iment n{tp k ) individual measurements are performed at phase tp k , k = 0,l, . . . , N — 1, 
then (F) can be estimated as 

N-l I n{tp k ) 

(F<«t>(JV)> = - £ -— - £ K F [x n (<p k ),<p k ], (189) 

iV fc=0 n Vfk) n= \ 

where x n (iph) is the result of the nth individual measurement at phase ip k - Taking the 
average over all estimates (F^^N)) yields (F(N)); i.e., the desired quantity within 

85 Equation (187) can be derived by substituting in eq. (186) for p(x, tp) the result of eq. (98) and 
recalling eqs. (92) and (100). Note that G£|™ = S m kS n i only holds for k — l~m — n. 

86 For the problem of measuring the state of single quantum systems, see Ueda and Kitagawa 
[1992]; Aharonov and Vaidman [1993]; Aharonov, Anandan and Vaidman [1993]; Imamoglu [1993]; 
Royer [1994, 1995]; Alter and Yamamoto [1995]; D'Ariano and Yuen [1996]. 
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the systematic error owing to phase discretization, 

N-l 

N 



(FM)(iV)) = (F(N)) = £ £ WW] (190) 

fc=0 



where 



v?jfc] = y da;p(a;, </? fe ) Xf(^, <Pk)- (191) 

Accordingly, the statistical fluctuation of (F^ cst \N)) in terms of the averaged esti- 
mates of the variance of (F^ est \N)) can be given by, on taking into account that the 
individual measuremenst are statistically independent of each other, 87 



7T 



4) = x E N-\tp k ){AK F [x(tp k ),tp k }}*, (192) 



k=i 
where 

{Ai^ F [x(v5 fc ), tp k ]} 2 = J dxp(x, ifk) K 2 F (x, ifk) - J dxp(x, tp k ) K F (x, (p k ) 



(193) 



is the variance of the sampling function, and N(<p k ) is the number of measurements 
per phase interval, N(tp k ) = n(<p k )N/ir. 

Let us mention that when the kernel K F (x, tp) is a strongly varying function of x in 
regions where p(x, tp) is non-negligible, then the first term in eq. (193) is much larger 
than the second one. In this case the second term can be neglected and the statistical 
error can be approximated by averaging the square of the kernel. Note that the same 
result is obtained if one assumes that the numbers of events yielding the values of 
x((pk) in given intervals of x (bins) are independent Poissonian variables (Leonhardt, 
Munroe, Kiss, Raymer and Richter, Th., [1996]; Leonhardt [1997c]; McAlister and 
Raymer [1997b]). However, these variables are neither strictly independent [their sum 
is always N(<p k )] nor Poissonian [probability that there is more than N(tp k ) events in 
a bin is zero]. Therefore, care must be taken before one decides to use the simplified 
error estimation. Whereas for the Fock-basis density matrix elements the simplified 
estimation is very good (the kernels are strongly oscillating), for the exponential 
moments of canonical phase one has to take into account both terms in eq. (193), 
otherwise the error would be overestimated (the kernels are slowly varying functions). 

From eq. (192) it can be expected that the statistical error depends sensitively on 
the analytical form of the sampling function. To give an example, let us consider the 
integral kernels K mn (x, tp), eq. (92), needed for sampling the density-matrix elements 
in the photon-number basis, Q mn . For chosen tp, K mn (x, tp) is an oscillating function of 
x, and with increasing distance d— \m — n\ from the diagonal the oscillations become 
faster and the oscillation range slowly increases. It turns out that with increasing m 
the diagonal-element variance (J 2 mm becomes independent of m; i.e., saturation of the 
statistical error for sufficiently large m is observed, o 2 Qmm <2/N (N, mean number 
of measurements per phase interval). On the contrary, the off-diagonal statistical 
error increases with d, without saturation. The influence of the quantum efficiency rj 
on the statistical error of the density-matrix elements is very strong in general. For 



87 



For simplicity, in eq. (192) it is assumed that F is a real quantity. 
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chosen m and n the oscillation range of K mn (x, ip; rj), eq. (96), increases very rapidly 
as 7] approaches the lower bound rj — 1/2, and the statistical error increases rapidly as 
well (for numerical examples, see D'Ariano, Macchiavello and Sterpi [1997]). 88 Using 
the experimentally sampled density-matrix elements for calculating the expectation 
values of other quantities, such as (F) according to eq. (160), the resulting statistical 
error is determined by the law of error propagation. It is worth noting that error 
propagation can lead to additional noise which is not observed if the quantities are 
also directly sampled from the measured data (provided that the sampling method 
applies). 

3.9.2. Least-squares method 

In the quantum-state reconstruction problems outlined in the foregoing sections a 
set of measurable quantities (in the following also referred to as data vector) is re- 
lated linearly to a set of quantities (state vector) that can be used to characterize 
the quantum state of the system under study. Both sets of quantities can be dis- 
crete or continuous or of mixed type. Typical examples are the relations (70) and 
(98), respectively, between the Wigner function and the density-matrix elements in 
the Fock basis, and the relations (112) and (140) between the Fourier components 
of the displaced Fock-state distributions and the atomic-state inversion, respectively, 
and the density-matrix elements in the Fock basis. A powerful method for inversion 
of such relations has been least-squares inversion. 89 The method has been used for 
quantum-state reconstruction in balanced optical homodyning (Tan [1997]), unbal- 
anced homodyning (Opatrny and Welsch [1997]; Opatrny, Welsch, Wallentowitz and 
Vogel, W., [1997]), cavity QED (Bardroff, Mayr, Schleich, Domokos, Brune, Raimond 
and Haroche [1996]; Bodendorf, Antesberger, Kim, M.S. and Walther [1998]) and for 
orbital electronic motion (Cline, van der Burgt, Westerveld and Risley [1994]). It 
has been used further to reconstruct the quantum state of the center-of-mass motion 
of trapped ions (Leibfried, Meekhof, King, Monroe, Itano and Wineland [1996]), the 
quantum state of a particle in an anharmonic potential and the quantum state of a 
particle that undergoes a damped motion in a harmonic potential (Opatrny, Welsch 
and Vogel, W., [1997c]). 

An advantage of least-squares inversion is that it is a linear method - the density 
matrix elements can be reconstructed in real time together with an estimation of the 
statistical error. Moreover, it allows for an easy incorporation in the reconstruction 
of various experimental peculiarities, such as non-unity quantum efficiency, finite 
resolution or discretization of the data, finite observation time, dissipative decay of 
the system, etc.. These aspects can hardly be treated on the basis of analytically 
determined (and existing) sampling functions. On the other hand, the method does 
not guarantee (similarly as any other linear method) that a reconstructed density 
matrix is exactly positive-definite (cf. § 3.9.3.). 

88 For the error of the exponential phase moments sampled from the quadrature-component dis- 
tributions in balanced homodyning, see Dakna, Opatrny and Welsch [1998], and for the error in 
quantum-state measurement via unbalanced homodyning and direct photocounting, see Wallen- 
towitz and Vogel, W., [1996a]; Banaszek and Wodkiewicz [1997b]; Opatrny, Welsch, Wallentowitz 
and Vogel, W., [1997]. 

89 The method of least squares was discovered by Legendre [1805] and Gauss [1809,1821] for solving 
the problem of reconstruction of orbits of planetoids from measured data. 
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To illustrate the method, let us assume that a distribution p(x, ip) of the type of 
a quadrature-component distribution is measured, and that p(x, p) can be given by 
a linear combination of all density-matrix elements g nn r of the quantum state to be 
reconstructed, with linearly independent coefficient functions S nn '(x,ip), 

P( X i V) = J2 S nn'(x, if) £W, (194) 

n, n' 

where S nn '(x, p>) need not be of the form used in eq. (98). 90 Since the density matrix 
of any physical state can be truncated at some value n max , the sum in eq. (194) is 
effectively finite. Direct application of least-squares inversion (Appendix D.) yields 
the reconstructed density-matrix elements q^) as 

Qn$ = [ dx [ dpK nn ,(x,p)p^ M \x,p), (195) 

where p^ M \x,p) is the experimentally measured distribution, X and $ being the 
intervals accessible to measurement. The integral kernel K nn i(x,p) is given by 

K nn i(x,<p)= F nn';mm'S^ m ,(x,ip), (196) 

m,m'<n max 

and F = G _1 , with the matrix G being defined by 

G mm >;nn> = dx dp S* mm , (x, (f) S nn > (x, (f) . (197) 
J X 

It can be proved by direct substitution that if the data correspond to the exact 
quantities p(x,ip), i.e., p( M \x,ip) —p(x,ip), then the reconstructed density matrix 
equals the correct one, i.e., q^J = Q nn >. On the other hand, if the experimental data 
suffer from some inaccuracies, then the reconstructed density matrix has the property 
that it reproduces the data as truly as possible (in the sense of least squares). In the 
above given formulas we have assumed that x and ip are continuos variables, and 
n is discrete. The formulas for other combinations of discrete and/or continuous 
arguments can be obtained in a quite similar way. 

An essential point of the method is the inversion of the matrix G, which requires 
the matrix to be sufficiently far from singularity; i.e., the data must carry enough in- 
formation about all the density-matrix elements that are desired to be reconstructed. 
Otherwise regularized inversion must be applied (Appendix D.). Regularized inversion 
usually decreases the statistical error of the reconstructed density-matrix elements, 
but on the other hand they are biased. Therefore, in practice such a degree of reg- 
ularization should be used for which the introduced bias is just below the statistical 
noise. 

3.9.3. Maximum-entropy principle 

As already mentioned, it is principally impossible to measure the exact expectation 
values of an infinite number of operators Aj in an expansion of the density operator 

90 In particular, when x is the position of a moving particle and ip corresponds to the time t, and 
the particle undergoes damping, then the quantum state evolves according to a master equation 
whose solution then determines S nn '(x,t). 
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of the type given in eq. (69), because any realistic experiment can only run for a finite 
time. So far, the exact formulas have been applied to the analysis of the incomplete 
measurements including an estimation of the error made. However the question may 
arise of how to obtain an optimum result of reconstruction of a quantum state when 
only a finite number of (Aj) has been measured in the experiment. An answer can 
be given using the Jaynes principle of maximum entropy 91 (Buzek, Adam, G., and 
Drobny [1996a,b]; Buzek, Drobny, Adam, G., Derka and Knight [1997]). 

Let us assume that the expectation values (Aj) of n quantities Aj, i — l,...n are 
experimentally determined. 92 The set of measured quantities can be regarded as a 
measure of the realized observation level. Certainly, there are a number of potential 
density operators g, Tr g— 1, which are compatible with the experimental results, i.e., 

r Fr[gA i ) = {A i ), i = l,...n. (198) 

Among them that density operator is chosen that maximizes the von Neumann 
entropy 93 

S[q] = -Tr(glng). (199) 

Introducing Lagrange multipliers, the resulting density operator g = Qs takes the 
familiar (grand-canonical ensemble) form 

Qs = Z- 1 exp(-p\ i A^, (200) 

where 



Z = Tr 



eX P ( - 2 



(201) 



and represents a partially reconstructed (estimated) density operator on the cho- 
sen observation level. Substituting in eq. (198) for g the density operator g s from 
eq. (200), a set of n nonlinear equations is obtained for the calculation of the n La- 
grange multipliers Aj from the measured expectation values (Aj). Any incomplete ob- 
servational level can be extended to a more complete observational level, in principle, 
by including additional observables in the scheme, which is usually associated with 
a decrease of the entropy. However, since rather involved calculations are required 
to be performed, the method has been studied for reconstructing the quantum state 
of a radiation-field mode on particular (not very high) observational levels (Buzek, 
Adam, G., and Drobny [1996a,b]; Buzek, Drobny, Adam, G., Derka and Knight 
[1997]) and/or low-dimensional systems, such as spin states (Buzek, Drobny, Adam, 
G., Derka and Knight [1997]). 



91 See Jaynes [1957a, b]. Probability distributions (or density operators) describe our stage of 
knowledge about physical systems. If we do not know anything, we usually assign uniform dis- 
tributions to the quantities (or a multiple of the unity operator to the density operators). If we 
have partial knowledge, we choose such distributions which are as broad as possible and still reflect 
our stage of knowledge. A suitable measure of the breadth is the entropy, one therefore seeks for 
such distributions (density operators) which maximize the entropy under the condition that known 
quantities are reproduced. 

92 Note that the determination of (Aj) already requires an infinite number of individual measure- 
ments which cannot be realized during a finite measurement time. 

93 See von Neumann [1932]. 
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Since the expectation values (A\) cannot be measured with infinite precision, only 
estimates can be inserted into eq. (198), and therefore it can happen that the solution 
does not exist; i.e., the non-precisely measured averages are not compatible with any 
density operator. To overcome this problem, the method can be combined with least- 
squares minimization (Wiedemann [1996] 94 ); i.e., the sum of squares of differences 

n „ 

C{{\}) = E {<^> (M) " Tr[es({*i})M } (202) 

as a function of the parameters Aj is tried to be minimized, (Aj)( M "> being the non- 
precisely measured averages. 

3.9.4. Bayesian inference 

The statistical fluctuations of the data are taken into account in the Bayesian in- 
ference scheme (Helstrom [1976], Holevo [1982], Jones, K.R.W., [1991, 1994]; Derka, 
Buzek, Adam, G., [1996]; Derka, Buzek, Adam, G., and Knight [1996]). Let us as- 
sume that the system under consideration is prepared in a pure state that belongs 
to a continuous manifold of states in a state space Q and consider a repeated iV-trial 
measurement of observables Aj, % — 1, . . . ,n, with eigenvalues Aj i . The determina- 
tion of the quantum state of the measured system is then performed in a repeated 
three-step procedure: 

(i) As a result of a (single) measurement of A,- L a conditional probabilty p(Aj.\g) 
is defined which specifies the result A,-, if the measured system is in state g = 

p(A ji \Q)=Tr(P A j), (203) 



where P A = \A jz )(A 



n 



(ii) A probability distribution po(g) defined on the space Q is specified such that 
it describes the a priori knowledge of the state to be reconstructed. The joint 
probability distribution p(Aj { , g) is then given by 

p(A j .,Q)=p(A ji \Q)p (Q). (204) 

When no initial information about the measured syste is available, then the 
prior probability distribution po(g) is chosen to be constant. 

(iii) Finally, the Bayes rule 95 is used to obtain the probability p(g\Aj.) of the system 
being in state g under the condition that Aj. is measured, 



94 In this paper, which is unfortunately unpublished, the operators Ai arc identified with the 
photon number n and quadrature-component projectors \xi, y>k){xu<fik\, where the subscript I labels 
a finite subset of the continuous quadrature-components at chosen phase ipt- Due to computational 
limits 4 phases and 13 values of x at each phase are considered. The partial reconstruction of the 
quantum state (in phase space) from computer-simulated homodyne data is performed for various 
states and yields results that reflect typical properties of the states sufficiently well. 

The Bayes rule expresses the relation between the conditional probabilities p(x\y) and p(y\x) as 
p(x\y)p(y) =p(x,y) =p(y\x)p(x), and hence p(y\x) =p(x,y)/p(x) =p(x,y)/ J dyp(x,y). 
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Now the procedure can be repeated, making a second measurement of some observable 
Ak (Ak = Ai or Ak ^ AA and using p(£|AjJ obtained from the first measurement as 
the prior probability distribution for the second measurement. Proceeding in the way 
outlined, the iV-trial conditional probability distribution p(g\{Aj.}) is given by 

n fx\fA \\ £ (g) Po(e) ( 0nf * 

M{A '' ]) = l^rmMf)' (206) 



where 



C(g) = fining) (207) 

i=i 



is the likelihood function (regarded as a function of g; the measured values {Aj t } 
playing the role of parameters). Note that C(g) is the probability distribution of 
finding the result {A,-.} in the sequence of N measurements under the condition that 
the state of the system is g. The partially reconstructed density operator g B is then 
taken as the average over all the possible states g, ge 

g B = J^dQp(g\{A n })g. (208) 

To apply the method, the state space f2 of the measured system and the cor- 
responding integration measure dO must be defined and the prior probability po(g) 
must be specified. In particular, the integration measure has to be invariant under 
unitary transformations in the space Q. This requirement uniquely determines the 
form of the measure. It should be pointed out that this is no longer valid when Q 
is tried to be extended to mixed states. Since a system which is in a mixed state 
can always be considered as a subsystem of a composite system that is in a pure 
state, the Bayesian reconstruction can be applied to the composite system and trac- 
ing the resulting density operator over the degrees of freedom of the other subsystem 
(Derka, Buzek and Adam, G., [1996], Derka, Buzek, Adam, G., and Knight [1996]). 
As already mentioned, the prior probability po(g) can be chosen constant if there 
is no a priori information about the state to be reconstructed. 97 It turns out that 
with increasing number of measurements the method becomes rather insensitive to 
the prior probability. In particular, when the number of measurements approaches 
infinity, iV — > oo, eq. (208) corresponds, on chosen observation level, to the principle 
of maximum entropy on a microcanonical ensemble. 98 Contrary to least-squares in- 
version, Bayesian inference always yields a (partially) reconstructed density operator 
which is (formally) positive definite and normalized to unity. However the price to 
pay is a rather involved procedure that has to be carried out in practice. For this 
reason, the method has been mostly considered - similar to the method of maximum 
entropy - for spin systems; i.e., for systems with small dimension of the state space 



96 Note that qb can correspond to a mixed state even though it is assumed that the system is 
prepared in a pure state. 

97 Note that when the integration measure dfi is not uniquely defined, then a prior probability 
that is constant with respect to a chosen integration measure need not be constant with respect to 
another one. 

98 Note that on a chosen (incomplete) observational level the two methods yield different fluctua- 
tions of the observables in general. 
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(Jones, K.R.W., [1991], Derka, Buzek and Adam, G., [1996], Derka, Buzek, Adam, 
G., and Knight [1996]). 

Another statistical method closely related to the Bayesian reconstruction method 
is the quantum state estimation based on the maximization of the likelihood function 
C(g), eq. (207). Here, that state g in the state space is selected for which C(g) attains 
its maximum. Because of the difficulty of finding the maximum of C(g) in higher- 
dimensional state spaces, a procedure was proposed which is based on a sequence of 
general inequalities satisfied by the likelihood function (Hradil [1997]). In this way, 
the problem can be tranformed to that of the diagonalization of an operator given by a 
linear combination of the projectors Paj. , where the expansion coefficients must solve 
a set of nonlinear algebraic equations. The method simplifies the Bayesian treatment 
but still guarantees the positive definiteness. Significantly, all the solutions based on 
the deterministic relation (69) between counted data (frequencies) and the desired 
density matrix are involved. Whenever such solution exists as a positively defined 
density matrix, then it should maximize the likelihood function as well. The method 
was applied successfully to the reconstruction of (low-dimensional) density matrices 
in the photon-number basis of a radiation-field mode from computer-simulated homo- 
dyne data, and a comparison with direct inversion of the linear basic relation between 
the measurable quantities and the density-matrix elements was given (Mogilevtsev, 
Hradil and Perina [1997]). 

4. Quantum states of matter systems 

In the preceding sections we have considered phase sensitive measurements of ra- 
diation fields and methods for reconstructing the quantum state of the fields from 
the measured data. The problem of quantum-state measurement and reconstruction 
has also been studied for various matter systems. Different matter systems require, 
in general, different detection schemes for measuring specific quantities that carry 
the full information on the quantum state of the system. Although these methods 
may be, at first glance, quite different from the methods outlined in § 2. for phase- 
sensitive measurements of light, there have been a number of analogies between the 
reconstruction concepts for radiation and matter. 

4.1. Molecular vibrations 

It was shown and demonstrated experimentally that the quantum state of molecu- 
lar vibrations can be determined using a tomographic method (Dunn, Walmsley and 
Mukamel [1995]) which resembles the one for a light mode outlined in § 3.1.. The 
method, called molecular emission tomography, is based on the fact that the time re- 
solved emission spectrum of a molecule allows one to visualise the time dependence of 
a vibrational wave packet within the excited electronic state from which the emission 
originates (Kowalczyk, Radzewicz, Mostowski and Walmsley [1990]). Alternatively, 
the desired information on the wave-packet dynamics can be obtained by photoelec- 
tron spectroscopy (Assion, Geisler, Helbing, Seyfried and Baumert [1996]). 

Let us assume that the molecule is prepared in a given vibrational quantum state 
in the excited electronic state. As can be seen from Fig. 17, for appropriately displaced 
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potential energy surfaces of the molecule the position of the vibrational wave packet 
can be effectively mapped onto the frequency of the emitted light. This fact is used 
for the tomographic reconstruction of the vibrational wave packet by measuring the 
time-resolved emission spectrum with a time resolution that is fast compared with 
the characteristic time period of the molecular state to be studied. The experimental 
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Figure 17: The vibronic energies for the A 1 !!^ —>■ transition of Na2 

clearly show the possibility to display the vibrational motion (in the excited 
state) in the time-resolved emission spectrum. (After Kowalczyk, Radzewicz, 
Mostowski and Walmsley [1990].) 

realization has been performed as follows (Dunn, Walmsley and Mukamel [1995]). A 
sample of Na2 molecules is illuminated by a 4 kHz train of laser pulses of 60 fs duration 
and mean wavelength of 630 nm. The laser pulses generate vibrational wave packets 
in the A*£+ state of the sodium dimer, which evolve with a time period of 310 fs. A 
fraction of the pulses is split off and plays the role of a time-gate shutter. The light 
emitted from the molecular sample is collected and focused synchronously with the 
split-off part of the exciting pulse onto a nonlinear crystal. A prism monochromator 
is used to filter the resulting sum-frequency and the field is recorded by a photon- 
counting photomultiplier. The resulting temporal resolution of the device is about 
65 fs. 

4.1.1. Harmonic regime 

The first reconstruction of the quantum state of molecular vibrations from a time- 
resolved emission spectrum was based on the assumption that only low vibrational 
quantum states are excited such that the relevant potentials can be approximated 
by harmonic ones. Furthermore, it was assumed that the vibrational frequencies in 
the two electronic quantum states, which contribute to the emission spectrum, are 
nearly equal. In this case, the vibronic coupling is caused solely by the displacement 
of the equilibrium positions of the potentials in the two electronic states. When these 
approximations are justified, then the time-gated spectrum 5(0, T) can be related to 
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the s-parametrized phase-space distribution P(q,p; s) = 2 1 P[a = 2 x ' 2 (g + ip); s] as 
(Dunn, Walmsley and Mukamel [1995]) 

S{Q,T) = J dyP[x(Q)cos{uT)+ysm{uT),ycos{uT)-x(n)sm{uT);s). (209) 

Here, Q and T are the setting frequency of the spectral filter and the setting time 
of the time gate, respectively (z/, vibrational frequency). The function x(Q) = (Q — 
k 2 u) I \^2ku describes the mapping of the position of the wave packet onto the emitted 
frequency, k being the ratio of the displacement of the potentials to the width of the 
vibrational ground state. The ordering parameter s in eq. (209) is related to the 
temporal width T _1 of the time gate as s = —{T/kv) 2 . 

Obviously, the relation (209) between the phase-space distribution P(q,p;s) of 
the molecular vibration and the measured spectrum S(Q, T) is very close to the basic 
relation (70) of optical homodyne tomography. Thus P(q,p; s) can be obtained from 
S(Q, T) by means of inverse Radon transform (see § 3.1.). A typical experimental 
result is shown in Fig. 18. 
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Figure 18: Tomographic reconstruction of an s-parametrized phase-space 
function for a vibrational wave packet in sodium from a measured set of 
emission spectra. The time gate duration was 65 fs, which implies s= —0.8. 
(After Dunn, Walmsley and Mukamel [1995].) 



4.1.2. Anharmonic vibrations 

In general, molecular vibrations are known to be significantly anharmonic when their 
excitations are not restricted to very small numbers of vibrational quanta. In such 
cases one cannot apply the approximate reconstruction procedure based on eq. (209). 
Due to the anharmonicity effects it is no longer possible to reconstruct the quantum 
state only from one half of a vibrational period. 

It was proposed (Shapiro, M., [1995]) to reconstruct the vibrational quantum 
state from the time-resolved (spectrally integrated) intensity I(T) = J dQ S(Q, T) of 
the light emitted by the molecular sample, which can be related to the density-matrix 
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elements g mn of the vibrational mode in the excited electronic state as 

I(T) -EfwE \ km \* kn exp(-iu mn T). (210) 

m,n k 

Here, v mn are the vibrational transition frequencies in the excited electronic state, 
and \ km ~ (^Ti) 2 1(^^)2 is determined by the Franck-Condon overlap i(/c|m) 2 of 
the vibrational wave functions in the two electronic states and the vibronic transition 
frequency u™* . Equation (210) reveals that, as long as the transition frequencies v mn 
are nondegenerate, the corresponding density-matrix elements g mn can be obtained, 
in principle, from an analysis of I{T) as a function of T. However, the separation of 
the density-matrix elements from each other may require a rather long time series. 
Further, the dimension of the set of equations to be inverted can be large, because of 
the large number of density-matrix elements that may contribute to the intensity of 
the emitted light. In the degenerate case, which is observed for the diagonal elements 
of the density matrix and for some off-diagonal elements due to the anharmonicities, 
additional information is needed. It was proposed to use the stationary spectrum of 
the light, whose determination requires an additional measurement. Alternatively, 
the two measurements can be combined such that the density-matrix elements are 
reconstructed from the time-resolved spectrum (Trippenbach and Band [1996])." 

It is worth noting that the dimension of the sets of equations that must be inverted 
numerically can be reduced substantially by employing the full information inherent 
in the time-resolved spectrum: 

S(n,T) ~ Yl Qrnn Y X kmKn exp(-iu mn T) g(Q ^"O (211) 

m,n k 

(Waxer, Walmsley and Vogel, W., [1997]). In eq. (211), the blurring function g{u) = 
exp(— u 2 /4T 2 ) is determined by the resolution time T _1 of the time gate (see § 4.1.1.). 
In practice, the time-dependent spectrum is available only in a finite time interval of 
size r, which can be taken into account by multiplying S(Q, T) by the corresponding 
sampling- window function G(T, r) in order to obtain S'(Q, T) = S(Q, T)G(T, r). The 
Fourier transform of S"(f2,T) with respect to T reads 

S'(Sl,v) ~ Y (212) 

m,m k 

with Giv — z/ m „,r) being the Fourier transform of the sampling window. Fixing 
the frequency v of the time-series spectrum, the structure of the window function 
ensures that only some of the density matrix elements contribute at this frequency, 
depending on the chosen size r of the sampling window. Let us assume that the 
density matrix elements g mini (i = 1, 2, 3, . . . , N) contribute to the spectrum at chosen 
frequency v. By choosing iV frequencies Q = Qi of the emission spectrum, one gets 
from eq. (212) a linear set of iV equations of rather low dimension that can be inverted 
numerically. 100 In this way, the method allows one to reconstruct rather complicated 

"Trippenbach and Band [1996] also discussed the inclusion of molecular rotations in the 
reconstruction. 

100 Note that the emission frequencies fti are determined by the two vibrational potentials involved 
in the vibronic emission, whereas the degenerate values of the vibrational frequencies v mn are deter- 
mined solely by the vibrational potential in the excited electronic state. The (within the resolution 
of the sampling window) degenerate density-matrix elements usually contribute to the emission 
spectrum at distinct frequencies f^. 
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quantum states on a time scale that can be shorter than the fractional revival time, 
as was demonstrated in a computer simulation of measurements (Waxer, Walmsley 
and Vogel, W., [1997]). Note that a reduction of the dimension of the set of equations 
to be inverted improves the robustness of the method with respect to the noise in the 
experimental data. 

It was also proposed to reconstruct the density matrix from the time-dependent 
position distribution according to eq. (102) (Leonhardt and Raymer [1996]; Richter, 
Th., and Wunsche [1996a,b]; Leonhardt [1997a]; Leonhardt and Schneider [1997]). 
This requires a scheme suitable for measuring either the position distribution of molec- 
ular vibrations or another set of quantities that can be mapped onto the position at 
different times. Note that this is hardly possible in molecular emission tomography 
in general. Finally it was proposed to use a wave-packet interference technique (Chen 
and Yeazell [1997]; § 4.5.1.) for reconstructing pure vibrational states in the excited 
electronic state (Leichtle, Schleich, Averbukh and Shapiro, M., [1998]). 

4.2. Trapped-atom motion 

Since the first observation of a single ion in a Paul trap 101 (Neuhauser, Hohenstatt, 
Toschek and Dehmelt [1980]) much progress has been achieved with respect to laser 
manipulation of the quantized motion of single atoms in trap potentials. Such sys- 
tems are of particular interest since the quantized low-frequency (~ MHz) motion is 
very stable, and laser manipulations allow one to prepare very interesting nonclassi- 
cal states. Until now, motional Fock states and squeezed states (Meekhof, Monroe, 
King, Itano and Wineland [1996]) as well as Schrodinger-cat type superposition states 
(Monroe, Meekhof, King and Wineland [1996]) have been realized experimentally. 

One might expect that the reconstruction of the quantum state of the center-of- 
mass motion of a trapped atom may be very similar to the reconstruction of molecular 
vibrations. However, the vibronic couplings in the two systems are basically differ- 
ent. In the case of a molecule, the vibrating atoms are close together within atomic 
dimensions and the change of the electronic state substantially alters the potential of 
the nuclear motion. In the case of an atom in a Paul trap the potential of the center- 
of-mass motion is externally given by the trap. In this case, electronic transitions can 
hardly affect the potential. The vibronic interaction in this system is induced by the 
interaction with radiation. Therefore one may expect that appropriate interactions 
of a trapped atom with laser fields may open various possibilities for measuring the 
motional quantum state. 

4.2.1. Quadrature measurement 

The first proposals to reconstruct the motional quantum state of a trapped atom 
were based on measuring the quadrature components of the atomic center-of-mass 
motion in the (harmonic) trap potential. In the scheme in Fig. 19 a weak electronic 
transition \g) = |1) <-> |e) = |2) of the atom is proposed to be simultaneously driven 
by two (classical) laser beams whose frequencies u r and u>b, respectively, are tuned 
to the first motional sidebands, u r = UJ2\ — v and Ub = 0021 + v of the electronic 
transition of frequency cu 2 i (Wallentowitz and Vogel, W., [1995,1996b]). Since the 



11 For the trap, see Paul, W., Osberghaus and Fischer [1958]. 
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linewidth of the transition is very small, the motional sidebands can be well resolved. 
For a long-living transition and in the Lamb-Dicke regime 102 these interactions are 
well described by Hamiltonians of the Jaynes-Cummings (and anti Jaynes-Cummings) 
type. 103 Assuming that the Rabi frequencies of the two (classically treated) laser fields 
are equal, the two Jaynes-Cummings interactions can be combined to an electron- 
vibration coupling that in the interaction picture reads as 

H' = ±h{l L V2(a- + a + )x(tp), (213) 

where Ql is the vibronic Rabi frequency of the driven transitions, and tp is the phase 
difference of the two lasers which can be controlled precisely. Assuming that the 
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Figure 19: Scheme of a trapped ion with a weak transition |1) <-> |2) and 
a strong transition |1) <-> |3). Two incident lasers of frequencies uo r = uj 2 i — 
v and Ub = o^2i + v are detuned from the electronic transition by the vibra- 
tional frequency v to the red and blue respectively. The laser driving the 
strong transition is used for testing the ground state occupation probabil- 
ity by means of resonance fluorescence. (After Wallentowitz and Vogel, W., 
[1995].) 




atom is initially prepared in the electronic ground state, P g (t)\ t =o = 1 (incoherent 
preparation), it can be shown that the atomic-state inversion at time t is nothing 
but the real part of the (scaled) characteristic function fy(x, tp) of the quadrature- 
component distribution p(x, tp) at phase tp of the center-of-mass motion, 

AP (mc) ^ = _ Re q ^y/2Q L t, tp) (214) 

102 In this regime, the (one-dimensional) spread of the motional wave packet, Ax, is small compared 
with the laser wavelength Al over 27r, Ax <C \l/2it. The Lamb-Dicke parameter, tild = 2tt(Ax)o/Xl, 
is a measure of the spread of the motional wave packet in the ground state of the trap potential 
relative to the wavelength. 

103 For the description of the dynamics in the Lamb-Dicke regime by a Jaynes-Cummings type 
Hamiltonian, see Blockley, Walls and Risken [1992]; Wineland, Bollinger, Itano, Moore, and Heinzen 
[1992], Cirac, Blatt, Parkins and Zoller [1994]. Note that the fact that the method is based on the 
Lamb-Dicke regime is not a serious limitation. In recent experiments (Meekhof, Monroe, King, 
Itano and Wineland [1996]) the frequency of the weak transition is in the GHz range and it is driven 
by two lasers in a Raman configuration. In this case the Lamb-Dicke parameter can be changed 
easily by the laser beam geometry. 
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[AP = P e — Pg\. Similarly, the imaginary part of the characteristic function can be 
obtained by preparing the atom initially in an appropriate coherent superposition of 
the two electronic quantum states such that 

V(y/2n L t, tp) = -AP (inc) (t) - AP (coh) (t). (215) 

Hence, varying the phase difference between the two driving laser beams, the 
complete information about the motional quantum state can be determined from the 
time evolution of the atomic-state inversion of the weak transition. This occupation 
probability can be probed by testing an auxiliary, strong transition |1) <-> |3) for the 
appearance of fluorescence (cf. Fig. 19). The existence and the absence of fluorescence 
detects the weak transition in state |1) and |2), respectively. This method has an al- 
most ideal quantum efficiency which is a great advantage for the present purpose. 104 
Applying eqs. (85) - (87), the reconstruction of density matrix of the center-of-mass 
motion in a quadrature-component basis from the measured signal can be accom- 
plished with a single Fourier integral, and application of eqs. (105) and (106) offers 
the possibility of direct sampling of the density matrix in the Fock basis from the flu- 
orescence signal. Further, using eqs. (B.24) and (B.27), s-parametrized phase-space 
functions can be inferred from the signal by performing a two-fold integral. 

The latter was demonstrated in another proposal suitable for measuring fy(z,tp) 
(D'Helon and Milburn [1996]). Here, the |1) <-> |2) transition (Fig. 19) is driven by a 
standing-wave laser pulse tuned to u 2 i whose duration Tp is much shorter than the 
vibrational period, ut p <C 1, and it is assumed that the center of the trap potential 
coincides with a node of the standing wave. In the interaction picture, the interaction 
Hamiltonian in the Lamb-Dicke regime then reads as 105 

H'=±m L (t)V2(a- + &+)x(t). (216) 

Here, x(t) corresponds to the freely evolving position operator of the center-of-mass 
motion, which for the harmonic motion considered is simply given by the quadrature- 
component operator at phase tp = ist, and the time- dependent (slowly varying) Rabi 
frequency flp{t) includes the shape of the pulse. Owing to the shortness of the driving 
pulse it can be assumed that its action at time t gives rise to a unitary kick described 
by the operator 



K = T exp 



-ih' 1 J dt' H\t') w exp \-i\6(&- + a+)x(t)] , (217) 



where — J dt' Q(t') is the pulse area. Obviously, its action on the state is formally 
the same as the action of the unitary operator that corresponds to the interaction 
Hamiltonian (213). Hence, A^ mc )p(t) and A( coh )p(t) measured immediately after 
application of the pulse at time t are again related to ^(z, tp) according to eq. (215), 
but now with 6 and vt in place of f^t and tp, respectively. 106 



104 For the experimental detection of a weak transition via fluorescence from another, strong tran- 
sition, see Nagourney, Sandberg and Dehmelt [1986]; Sauter, Neuhauser, Blatt and Toschek [1986]; 
Bergquist, Hulct, Itano and Wineland [1986]. 

105 Notc that the Hamiltonian (216) is not based on a vibrational rotating wave approximation as it 
is the case for the Hamiltonian (213). In eq. (216) the position operator results from the Lamb-Dicke 
approximation of the (operator-valued) mode function of the standing wave, sin(fc^i) w k^x. 

106 Whereas in the scheme of Wallentowitz and Vogel, W., [1995] the phase difference between the 
two laser beams must be controlled, the scheme of D'Helon and Milburn [1996] requires control of 
short-pulse areas. 
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It is well known that a squeezed coherent state approaches a quadrature-compo- 
nent eigenstate as the strength of squeezing goes to infinity; i.e., 

\x,<p)cx \im D(xe ilp /V2)S{\t\e i2v )\0), (218) 

|£|-»oo 

with D(ot) the coherent displacement operator, and S(£) = exp[(£*a 2 — £a^ 2 )/2] the 
squeeze operator. 107 Hence, the quadrature-component distribution p(x,<p) can be 
asymptotically given by 

p(x,<p)<x Jim (0|£'|0), (219) 

where g' is a coherently displaced and squeezed version of the density operator g to 
be determined, 

g> = 5 t (|^|e i2v ) D^xe^/VS) gD(xe ilp /V2) S(\Z\e i2 v). (220) 

Equations (219) and (220) reveal that when the quantum state to be measured can be 
squeezed strongly and displaced coherently, then the quadrature-component distribu- 
tions can be obtained, in principle, from the occupation probability of the motional 
ground state that is measured after these manipulations, the phase being controlled 
by the free evolution of the vibrating system (Poyatos, Walser, Cirac, Zoller and Blatt 
[1996]). To realize the scheme, the following procedure was proposed: 

(i) Wait for a time t such that ip = ut. 

(ii) Perform a sudden displacement of the center of the trap to the right for a 
distance d such that x = ^mv/ (2h) d, m being the mass of the atom. 

(iii) Change the trap frequency instantaneously from v to (lower) z/'. 108 

(iv) Determine the population of the motional ground state. 109 

Knowing p(x, ip) for all phases p in a n interval, the complete information on the 
motional quantum state is available and the (tomographic) methods of quantum- 
state reconstruction outlined in § 3. applies. 110 



4.2.2. Measurement of the Jaynes— Cummings dynamics 

Let us consider the dynamics of a trapped ion that is driven by a single (classical) laser 
beam in a resolved sideband regime 111 and assume that the laser frequency is tuned 
to the (red) /cth-order motional sideband of the weak electronic transition |1) <-> |2), 

107 For a proof, see, e.g., Vogel, W., and Welsch [1994]. 

108 The sudden change of the trap frequency leads to squeezing with the squeeze parameter being 
|£| = ln(v/v')/2 (Janszky and Yushin [1986]). Note that a significant change of the trap frequency 
is needed since the measurement scheme requires strong squeezing. 

109 It was demonstrated experimentally that the motional number statistics can be determined 
from the measured Jaynes-Cummings revivals (Meekhof, Monroe, King, Itano and Wineland [1996]; 
§ 4.2.2.). 

110 The steps (i), (ii) and (iv) were proposed in order to measure the Q function. 
m In the practical realization the single laser may also be replaced with two lasers driving a dipolc- 
forbidden transition in a Raman configuration, which essentially yields the same basic Hamiltonian. 
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W21 — ku. In this case the resulting coupling (in the interaction picture) between 
the electronic transition and the center-of-mass motion has the form of a nonlinear 
multiquantum Jaynes-Cummings interaction (Vogel, W., and de Matos Filho [1995]), 

H' = ±htt L a + f ik) {a^a) a k + H.c, (221) 
where the operator function reads in normally ordered form as 

00 (vr, \2l + k 

/^(^^e-^E^y^^^ (222) 

(t)ld, Lamb-Dicke parameter). It describes the influence of the interference between 
the driving laser wave and the (extended) atomic wave function on the dynamics 
of the system. The nonlinear Jaynes-Cummings model can be solved exactly. In 
particular, when the atom is initially prepared in the excited electronic state, then 
the atomic-state inversion is given by an equation of the type of eq. (137), 112 

00 

AP(t) = £ g nn cos(Q nn+k t) , (223) 

n=0 

where now Vt nn+ k = Vti J {n\f <yk ^a k \n+ k), and g nn are the motional density-matrix ele- 
ments at the initial time t — 0. 

The motional quantum state can be reconstructed from the measured (nonlinear) 
Jaynes-Cummings dynamics by applying typical methods outlined in § 3.. When the 
motional quantum state which is desired to be reconstructed is coherently shifted [g— > 
D' (a) qD (a)] before AP(t) is measured, then the Fock-state probability distribution 
Pn = Qnn in eq. (223) must be replaced with the displaced Fock-state probability distri- 
bution p n (a) = g nn {a) = (n,a\g\n,a) [cf. eq. (138)], and inversion of eq. (223) yields 
p n (a). Changing a in a succession of (ensemble) measurements, the complete infor- 
mation about the original (unshifted) motional quantum state can be inferred from 
the measured data. In particular, applying the method outlined in § 3.5., pointwise re- 
construction of phase-space functions is feasible. This was demonstrated successfully 
by reconstructing the Wigner function [eq. (130)] of a single 9 Be + ion that is stored 
in a rf Paul trap (Leibfried, Meekhof, King, Monroe, Itano and Wineland [1996]). 
In the experiment, the relevant oscillation frequency in the trap potential is z//27r~ 
11.2 MHz, and the transition between the states \g) and |e) is a stimulated Raman 
transition between the hyperfine ground states 2 Si/ 2 (F = 2, m F = —2) and 2 Si/ 2 (F 
= 1, m F = —1), respectively, which are separated by about 1.25 GHz. The coherent 
displacement of the initially prepared motional quantum state is realized by applying 
a classical, spatially uniform rf field. In Fig. 20 an example of a reconstructed Wigner 
function of a motional number state is shown. 

The density matrix in the Fock basis, g mn , can be obtained from the displaced 
Fock-state probability on a circle, applying the method outlined in § 3.3.2.. This 
was also demonstrated in the above mentioned experiment (Leibfried, Meekhof, King, 
Monroe, Itano and Wineland [1996]). In Fig. 21 an example of a reconstructed density 
matrix of a superposition of two Fock states is shown. 



112 The dynamics were realized experimentally and used to determine the (initial) excitation statis- 
tics g nn of the motional state by inverting eq. (223) (Meekhof, Monroe, King, Itano and Wineland 
[1996]). The same dynamics can be realized by tuning the laser onto the corresponding blue sideband, 
W21 =cji + kv, and preparing the electronic subsystem initially in the ground state. 
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Figure 20: Surface and contour plots of the reconstructed Wigner function 
W(a) of the |^ = 1) motional number state. The plotted points are the result 
of fitting a linear interpolation between the actual data points to a 0.1 by 
0.1 grid. The octagonal shape is an artifact of the eight measured phases 
per radius. (After Leibfried, Meekhof, King, Monroe, Itano and Wineland 
[1996].) 



The method of coherently displacing the motional quantum state to be detected 
can be extended in order to measure the Wigner function rather than reconstructing 
it from the measured data (Lutterbach and Davidovich [1997]). 113 For this purpose 
the laser driving the Jaynes-Cummings dynamics is tuned on resonance to the elec- 
tronic transition, so that the Hamiltonian (221) for k — applies. For sufficiently 
small Lamb-Dicke parameter, <C 1, the nonlinear function in eq. (222) can 
be approximated by ~ 1 — r]\ D (a)a + 1/2). Consequently, the vibronic Rabi fre- 
quencies read as il nn Q, l [1 — rf LD (n + 1/2)]. Inserting this result into eq. (223), 
with Q nn {a) in place of g nn , and choosing the interaction time as t — tw — ^ / (^Vld) 
yields 

oo 

AP(t w ) » sin(7r/77l D ) ]T g nn (-l) n = ^sm(n/rjl D ) W(a), (224) 

71=0 

which reveals that under the conditions given the measured electronic-state occupa- 
tion probability directly gives the Wigner function at the phase-space point a. 

Further it was proposed to reconstruct the density matrix in the Fock basis follow- 
ing a line similar to that given by eqs. (134) - (136) for cavity QED (Bardroff, Leichtle, 
Schrade and Schleich [1996]). In cavity QED the /c-photon transitions needed to cal- 
culate the off-diagonal density matrix elements g nn +k from the measured atomic-state 
inversion (provided that initially a coherent superposition of the two states is pre- 
pared) can hardly be realized for larger values of k. In the case of a trapped ion the 
situation is improved significantly. The /c-quantum interaction of Jaynes-Cummings 
type as given by eq. (221) can be realized as outlined above, by tuning the laser on 



113 For comparison, see the approach to direct measurement of the Wigner function in cavity QED 
as outlined in § 2.4.. 
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Figure 21: Reconstructed density-matrix amplitudes \p\ = \g\ mn of an approx- 
imate l/y/2(\n — 0) —i\n — 2)) motional state. The state was displaced by 
\a\ = 0.79 for eight phases on the circle. (After Leibfried, Meekhof, King, 
Monroe, Itano and Wineland [1996].) 



resonance to the kth motional sideband. This enables one to reconstruct the density 
matrix elements g nn +k f° r all k in principle, by successively tuning the laser on all 
sidebands. Needless to say, this requires sufficiently large values of the Lamb-Dicke 
parameter, since otherwise the coupling of the laser to high-order sidebands is very 
small. Formally, the difference to cavity QED only consists in the frequencies Q nn+ k, 
which are different for the two systems. 114 

It was also proposed to reconstruct the density matrix in the coherent-state ba- 
sis by combining coherent displacements with motional ground state measurements 
(Freyberger [1997]). The scheme allows a pure motional quantum state to be recon- 
structed by coherent displacements on a circle. Reconstruction of the density matrix 
of a mixed state requires additional variation of the displacement amplitude. A first 
displacement is introduced before a filtering measurement is performed in order to de- 
cide whether or not the atom is in the motional ground state. A second displacement 
is performed which is followed by detecting the probability that the atom is in the 
motional ground state, which may be obtained from the Jaynes-Cummings dynamics 
as outlined above. 

4.2.3. Entangled vibronic states 

Let us consider the situation where the motional state is entangled with the elec- 
tronic state of the system and ask for measuring the complete vibronic state. Such a 
measurement may be realized by the interaction Hamiltonian (221) for k = (Wallen- 
towitz, de Matos Filho and Vogel, W., [1997]). In the scheme the initially prepared 
vibronic state g is first displaced coherently in the phase space of the motional sub- 

114 At this point it should be noted that in a Paul trap the applied rf field may modulate the motion 
of a trapped atom, which gives rise to the so-called micromotion. The rf frequency is usually large 
compared to the (effective) motional frequency, and its effect may consist in a modification of the 
excitation-dependent function in the interaction Hamiltonian (221) (Bardroff, Leichtle, Schrade 
and Schleich [1996]). Apart from such modifications, the structure of the Hamiltonian is preserved. 
For details of the treatment of the micromotion, see the original article. 
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system, g — > g{a) = D\a)gD{a). Next the driving laser causing the electron-motional 
coupling is switched on for an interaction time r. This procedure is followed by prob- 
ing the atom for fluorescence on the strong, auxiliary transition. Provided the atom 
is detected in the excited electronic state |e) (no fluorescence), the density operator 
g(r) of the system reduces to 

&(T)~\e)(e\®Q M (T), (225) 

where £>m( t ) ~ ( e lf?( r )l e ) is the corresponding density operator of the motional sub- 
system. Its diagonal elements in the Fock basis read 

[QM{r)\ nn ~ Im[^"(a)] sin(fi n „r) 

+ g n e :(a) cos 2 (|fi n „r) + g™(a) sin 2 (ifi„ n r) . (226) 

Equation (226) reveals that appropriately chosen interaction times r allow the dis- 
placed vibronic density-matrix elements g^^a) to be mapped onto the reduced mo- 
tional number statistics [gM(j)} nn . It is worth noting that the underlying vibronic 
coupling is not only suited for this mapping but also for a QND measurement of 
the reduced motional number statistics, since the Hamiltonian commutes with the 
motional number operator. 115 Note that the reduced motional number statistics can 
also be determined using other methods, such as the methods outlined in § 4.2.2.. 

Knowing the displaced vibronic density-matrix elements Q™(ot) as functions of a, 
the density matrix g™ b n of the original vibronic quantum state can be reconstructed, 
following the lines given in § 3.. In particular, the elements Q^(ct) can be summed up 
to obtain the vibronic quantum state in terms of the Wigner-function matrix W a b(a), 

r\ CO 

^(«) = -£(-l)>r(«) (227) 

(Wallentowitz, de Matos Filho and Vogel, W., [1997]). The Wigner-function matrix 
W ah (a) has the following properties. Its trace with respect to the electronic subsystem 
is the (reduced) motional Wigner function. Integrating W a b(a) with respect to the 
motional phase-space amplitude a yields the (reduced) electronic density matrix. A 
typical quantum state for which such a measurement scheme would be of interest is 
an entangled state of the type \ip) ~ ( 1 2) |ct) ± |1)| — a)). 116 Clearly, the nonclassical 
features of such a state are completely lost when one measures only the reduced 
motional quantum state. 



4.3. Bose— Einstein condensates 

The recent progress in evaporative cooling of an atomic gas has rendered it possible 
to realize Bose-Einstein condensation experimentally (Anderson, Ensher, Matthews, 

115 For the application of this type of vibronic coupling to QND measurement of the motional 
excitation of a trapped atom, see de Matos Filho and Vogel, W., [1996]; Davidovich, Orszag and 
Zagury [1996]. 

116 States of this type have been realized experimentally (Monroe, Meekhof, King and Wineland 
[1996]). For a reconstruction of W a b(a) from a computer simulation of measurements, see Wallen- 
towitz, de Matos Filho and Vogel, W., [1997]). 
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Wieman and Cornell [1995]; Bradley, Sackett, Tollett and Hulet [1995]; Davies, 
Mewes, Andrews, van Druten, Durfee, Kurn and Ketterle [1995]). The Bose-Einstein 
condensate (BEC) represents a macroscopic occupation of the ground state of the 
gas, which is an important signature of quantum-statistical mechanics. Moreover, 
coherence effects, such as interference fringes between two condensates, could be 
demonstrated (Mewes, Andrews, Kurn, Durfee, Townsend and Ketterle [1997], An- 
drews, Townsend, Miesner, Durfee, Kurn and Ketterle [1997]; for a theoretical in- 
terpretation, see Wallis, Rohrl, Naraschewski and Schenzle [1997]). In view of these 
feasibilities it is interesting to get more insight in the exact nature of the quantum 
state of BEC. 

Let us consider a two- mode BEC consisting of N atoms, 117 where the two modes 
may correspond to different hyperfine states. To reconstruct the quantum state, it 
was proposed to introduce a controllable phase shift (f> between the modes and to mix 
them by applying a (lossless) beam-splitter transformation. The numbers of atoms 
in the two modes are counted, and the two-mode quantum state is inferred from 
joint counting statistics (Bolda, Tan and Walls [1997], Mancini and Tombesi [1997b], 
Walser [1997]). 118 Neglecting collisions between the atoms, the two-mode density 
operator g is transformed as 



# = fr{P,<i>)QU(p,<l>), (228) 



with 



U(P, 0) = exp U/3 (a\a 2 e i(t> + 4&ie~^) . (229) 



[cf. eq. (7)], where the mode subscripts 1 and 2 label the corresponding hyperfine 
states. Raman transitions between these states can be made by optical pulses that 
are off-resonance from an excited state 3. An rf field may be used to couple state 1 to 
another state 4, and atoms in state 4 are repelled from the trap. The beam splitter 
transformation can be realized by applying two copropagating light pulses that are 
detuned from the 1—3 and 2 — 3 transitions. The phase is controlled by the phase 
difference of the lasers. The pulses are followed by an rf 7r-pulse that transfers the 
state 1 to 4. Consequently, the atoms in state 4 fall freely and those in state 2 remain 
trapped. Eventually, the seperate groups of trapped and untrapped atoms are counted 
(for the experimental setup, see also Mewes, Andrews, Kurn, Durfee, Townsend and 
Ketterle [1997], Andrews, Townsend, Miesner, Durfee, Kurn and Ketterle [1997]). 
The probability P m {4>) = -Pn,jv-n(0) of counting n and N — n atoms, respectively, in 
the modes 1 and 2 for a phase shift setting of <fi can be given by 119 

P m (<f>) = (m\U\(3A)QU{PA)\m) = E^W, (230) 

i,v 

where T^'(0) is definded by the corresponding matrix elements of the beam splitter 
transformation in the basis of the states \m). 



117 The two-mode BEC is considered in view of the lack of a coherent reference state within a 
condensate of a fixed number of atoms. 

118 The three proposals are basically similar. Concerning the practical realization of such a mea- 
surement, we briefly outline the rather detailed scheme proposed by Bolda, Tan, and Walls [1997]. 
For the interaction with radiation of an atomic Bose gas in an isotropic harmonic oscillator potential, 
see also Javanainen [1994]. 

119 For simplicity, the abbreviating notation |m) for the two-mode states |n)i ® \N — n) 2 is used. 
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The density matrix of the two-mode state can be obtained by numerically invert- 
ing eq. (230) (Bolda, Tan and Walls [1997], Mancini and Tombesi [1997b]) by means 
of the least-squares method (Appendix D.; see also § 3.9.2.). A further simplification 
can be achieved by taking the Fourier transform of P m {4>) with respect to <fi (Bolda, 
Tan and Walls [1997]). Since the sth Fourier component P^ = {2n)- 1 J 27r &<p e is ^ P m (<p) 
is related only to the matrix elements whose row and column indices differ by s, the 
corresponding blocks of equations can be inverted separately (cf. § 3.3.2.). Alter- 
natively, the inversion can also be performed analytically by applying orthogonality 
relations for Clebsch-Gordan coefficients (Walser [1997]). 

The complex amplitude of a BEC field, ip(r,t) = (ip(r,t)), is frequently assumed 
to satisfy a nonlinear Schrodinger equation 



(m, particle mass; U, potential; g, atom-atom interaction constant). The density 
of particles \ip\ 2 can be measured by phase-contrast imaging (Andrews, Mewes, van 
Druten, Durfee, Kurn and Ketterle [1996]). In order to determine ip(r,t) for given 
\i(j(r,t)\ 2 , it was proposed to use a nonlinear propagation for a trial function tp, so 
that the correct function ip can be found by adaptive modification of the trial function 
(Leonhardt and Bardroff [1997]). 

4.4. Atomic matter waves 

Stimulated by the successful experimental demonstration of atom interferometry 
(Carnal and Mlynek [1991], Keith, Ekstrom, Turchette and Pritchard [1991]) the 
study of coherence properties of atomic matter waves has been of increasing inter- 
est. In particular, several methods have been considered for determining the (single- 
particle) quantum state of both the transverse and the longitudinal center-of-mass 
motion of the atoms, and first experiments have been performed. 

4.4.1. Transverse motion 

The density matrix of the (with respect to the direction of propagation) transverse 
motion of a single particle can be determined by methods of refractive (atom) op- 
tics (Raymer, Beck and McAlister [1994], Janicke and Wilkens [1995]). Combining 
the effect of lenses and propagation, it is possible to reconstruct the quantum state 
from measured position distributions by tomographic methods (§ 3.1.). The desired 
rotation of the quadrature components can be realized by appropriately changing the 
parameters of the experimental set-up, such as the focal lengths of the lenses, their 
positions and the positions of the detectors (Raymer, Beck and McAlister [1994]). 120 
In the first experimental reconstruction of the transverse quantum state of an 
atomic matter wave (Kurtsiefer, Pfau and Mlynek [1997], Pfau and Kurtsiefer [1997]) 
the scheme is simplified by avoiding the use of lenses. Metastable helium atoms are 
used whose quantum state of the transverse motion is prepared by a combination of 




(231) 



i20p or pQggibie realizations of optical elements in atom optics, see, e.g., Adams, Sigel and Mlynek 
[1994]. 
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Figure 22: Measured time-resolved spatial atomic distribution in the Fraun- 
hofer diffraction limit. The two wave packets emerging from the individual 
slits overlap almost completely, and a maximal visibility of the interference 
fringes is observed. (After Pfau and Kurtsiefer [1997].) 



an entrance slit and a double slit, which realizes a nonclassical two-peaked atom- wave 
state. The atoms then propagate for some distance to a time- and space-resolving 
detector, which allows one to measure the transverse atomic distribution of the atoms 
with a spatial resolution down to 500 nm, and the arrival time of the atoms with an 
accuracy of 100 ns. This renders it possible to record spatial atomic distributions for 
different free-evolution times of the atom- wave packet after the double slit (Fig. 22), 
which are suited for reconstructing the quantum state. The measured distributions 
are position distributions of the sheared Wigner function 



(to, mass of the atom) of the freely evolving atom- wave packet after the double slit. 
Knowing the position distributions for a sufficiently large set of evolution times, the 
original quantum state can be obtained (§ 3.3.1.; for reconstruction of an original 
Wigner function from marginals of the sheared Wigner function, see also Lohmann 
[1993]). In particular, the measured position distributions correspond to the quadrature- 
component distributions whose phase parameter is given by 



(x , scaling length). Obviously, the phase ip can only be scanned between and n/2, 
and in a real experiment this interval is further reduced in general. Since a precise 
reconstruction of the quantum state requires a phase interval of size n, some a priori 
information about the state should be available (cf. § 3.9.1.). In the present case 
it may be possible to assume symmetry in position space for the Wigner function, 
because of the preparation of the quantum state by a double slit. The Wigner function 




(232) 
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reconstructed (with a limited range of phases if) from the measured diffraction pattern 
in Fig. 22 by applying inverse Radon transform (§ 3.1.) is shown in Fig. 23. 
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Figure 23: Wigner function of the atomic motional state immediately behind 
the double slit, reconstructed from the measured time- resolved spatial atomic 
distribution in Fig. 22. (After Pfau and Kurtsiefer [1997].) 



Although the shape of the Wigner function and especially the oscillating interfer- 
ence part are well reproduced, there is a systematic error owing to the limited range of 
phases accessible in the experiment. For example, the reconstructed Wigner function 
appears sheared. This shear could be avoided by forced symmetrization of the Wigner 
function in position. Using the Wittaker- Shannon sampling theorem (see, e.g., Marks 
[1991]), it can be proved which features of the reconstructed Wigner function of an 
atom-wave packet are obtained correctly when the set of data is not tomographically 
complete (Raymer [1997]). The analysis shows that the density matrix in the mo- 
mentum basis is correct everywhere except in an excluded region of width (<5p) m i n 
around the diagonal, with (5p) m i n being inversely proportional to the time interval 
over which the diffraction pattern is recorded. Transferring this result to the recon- 
structed Wigner function, it turns out that there is a low-frequency error in W(x,p) 
as a function of x, whereas the high-frequency behaviour is essentially correct. 

In another tomographic scheme for reconstructing the Wigner function of an ini- 
tially prepared transverse motional state, a set-up is considered, in which a classical 
standing light field that is strongly detuned from the atomic transition serves as a 
thick gradient-index lens for two-level atoms which are prepared in the ground state 
and cross the light field close to a node (Kienle, Fischer, Schleich, Yakovlev and 
Freyberger [1997]; for a set-up of this type, see also Fig. 8). It was found that here 
the phase if of the quadrature-component distribution associated with the measured 
time-resolved position distribution can be tuned over a complete ir interval. 

In a modified (interferometric) scheme it is assumed that the phase <fi of the mode 
function of the standing-wave laser field can be controlled via a movable mirror, which 
is used to create the standing wave (Freyberger, Kienle and Yakovlev [1997]; Kienle, 
Fischer, Schleich, Yakovlev and Freyberger [1997]). It is shown that when a pure 
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motional state is realized, then the phase of the free-evolving wave function can be 
inferred from the atomic position distributions for phases of the laser field y?L = — 7r/4, 
0,7r/4 and n/2 at chosen time. Knowing from a pre-measurement the absolute value 
of the wave function for this time, the initial wave function can then be reconstructed 
by inverting the free time evolution. 

4.4.2. Longitudinal motion 

In order to determine the quantum state of the longitudinal center-of-mass motion of 
atoms in an atomic beam, it was proposed to reconstruct the Wigner function tomo- 
graphically via state-selective time-dependent measurement of longitudinal position 
distributions (Kokorowski and Pritchard [1997]). In the scheme, two-level atoms that 
are initially prepared in the ground state are considered. Located at an adjustable 
position in the apparatus is an electromagnetic excitation region that is assumed 
to be effectively a n pulse; i.e., any particle exposed to the radiation in the region 
is completely excited. Downstream, a state-selective detector counts the number of 
particles which have made a transition into the excited state owing to the 7r pulse. 
Time dependence may be introduced either by operating the excitation region in short 
pulses at definite times, or by using an excited state whose lifetime is much shorter 
than the desired time resolution. In the latter case, detection consists of observing 
spontaneously emitted radiation, and the time is simply the time of detection. Mea- 
surement of the longitudinal position distribution as a function of time then yields, 
similar to eq. (233), the quadrature-component distributions for phase parameters in 
a 7r/2 interval. Again, symmetry assumptions may be made in order to compensate 
for the lack of accessible phases. 

It was also proposed to reconstruct the quantum state in terms of the density 
matrix in the energy basis by using an interferometric method that is based on a 
generalization of Ramsey's classic separated oscillatory fields technique 121 and the 
observation that the longitudinal momentum of an atom can be shifted coherently 
via interaction with off-resonant radiation (Dhirani, Kokorowski, Rubenstein, Ham- 
mond, Rohwedder, Smith, Roberts and Pritchard [1997]). In the interferometer an 
atom-wave in the internal (electronic) ground state is split into ground-state and 
excited-state components by an electromagnetic field (detuning Si) at position z\. 
The excited-state component receives a momentum shift h/Ski (~<5i/t> if the kinetic 
energy mv 2 /2 is much larger than hSi). The remaining ground-state component is 
split again by a second field (detuning S 2 ) at x 2 , where the excited-state component 
is momentum-shifted by hAk 2 . Both field regions are assumed to have different de- 
tunings (Si ^S 2 ). The excited state detected at the interferometer output contains a 
coherent superposition of two distinct energy components, each with a different corre- 
sponding longitudinal momentum. It can be shown that the probability of detecting 
an atom in the excited state at location z (output port) and time t is given by 122 



121 For details, see Ramsey [1956]. 

122 Here it is assumed that CI' — CI" <C CI' , CI"; for the relation without this restriction, see Dhirani, 
Kokorowski, Rubenstein, Hammond, Rohwedder, Smith, Roberts and Pritchard [1997]. 




X cos 



2 {\ [(Si - 5 2 )t - Ak[(z - zi) + Ak' 2 (z - z 2 )]} , (234) 
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where g(fl', fl") is the motional density matrix in the energy basis (strictly speaking, 
n = E/h). Equation (234) reveals that the measured excited-state probability P e (z,t) 
is determined by a convolution of the density-matrix elements g(fl', fl") with two 
factors. The exponential represents the free evolution of the density matrix, and the 
cosine squared contains the phase difference accumulated by the two paths of the 
interferometer. Eventually, eq. (234) can be inverted in order to obtain g(Q', Q") 
in terms of P e (z,t) by Fourier transforming twice after appropriate change of the 
variables (for details, see Dhirani, Kokorowski, Rubenstein, Hammond, Rohwedder, 
Smith, Roberts, and Pritchard [1997]). 

4.5. Electron motion 

So far knowledge of electronic quantum states has been typically required in order to 
reconstruct motional and vibrational quantum states of atoms and molecules. This 
may have brought the reconstruction of the motional quantum state of electrons into 
question (for the spin, see § 4.6.). 

4.5.1. Electronic Rydberg wave packets 

Highly excited Rydberg states have offered novel possibilities of preparing electron- 
wave packets. Since the early experiments at the end of the 1980s (Yeazell and Stroud, 
Jr., [1988]; ten Wolde, Noordam, Lagendijk and van Linden van den Heuvell [1988]) 
various schemes have been used in order to produce Rydberg wave packets (for a 
brief review of generation and detection of Rydberg wave packets, see Noordam and 
Jones, R.R., [1997]). A Rydberg wave packet is produced whenever several Rydberg 
states are excited coherently. In particular, it has also been possible to prepare an 
atomic electron in a coherent superposition of two wave packets which are localized 
spatially at the opposite extremes of a Kepler orbit (Noel and Stroud, Jr. [1996]). 
On a sufficiently short time scale these nonclassical states are close to the even/odd 
coherent states of the harmonic oscillator. To demonstrate important features of these 
states, two kinds of measurements have been used. State selective ionization allows 
one to verify that only every other atomic level is populated, and a Ramsey fringe 
measurement verifies the coherence of the superposition. 

A wave-packet interference technique was proposed to measure engineered atomic 
Rydberg wave functions (Chen and Yeazell [1997]). 123 It is assumed that the Rydberg 
wave packet is produced by (a sequence of) short optical (Gaussian) pulses whose 
interaction with the atomic system can be treated in lowest-order perturbation theory. 
Describing the atomic system by the radial wave function, 

*(r, t) = agMe-^-Xir) + £ a n (t)e-^ l u n {r), (235) 

n 

where a n and uo n are the Schrodinger amplitudes and eigenenergies for the Rydberg 
states and a g and u g are these for the ground state, the excited-state amplitudes after 
the last excitation pulse, a n , are given by 

a n = -±in n e- A * a2 / 2 z n , (236) 



123 The method was also proposed for reconstructing pure vibrational states in a molecule (Leichtle, 
Schleich, Averbukh and Shapiro, M., [1998]). 
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where z n = l + Yl i exp[i(u + A n )rj] (Q n , Rabi frequency of the excitation between the 
ground state and a given Rydberg state n; A n = uj n — u, detuning; a, pulse width; 
Tj, time delay for the ith individual pulse). When a probe pulse (identical to the 
previous pulses) interacts with the atom, then the new excited-state amplitudes 

a' n = -\Xl n e~ A ^ 2 / 2 {z n + exp[*(cu + A n )r]} (237) 

are produced (r, delay time for the probe pulse). By measuring the population 
in each eigenstate, \a' n \ 2 , for different phases between excitation and probe pulses 
(on applying high- resolution technique, such as state-selective field ionization), the 
complex quantities z n can be obtained; i.e., the (radial) wave function created by the 
excitation pulse may be inferred. Note that the pump and probe pulses need not be 
identical, but they must be phase coherent with each other. 

Using optical pulses, Rydberg wave packets can usually be detected when they 
are near to an ionic core. This limitation can be overcome by using subpicosecond, 
unipolar electromagnetic field (half-cycle) pulses, 124 which have been a powerful tool 
for studying dynamics in weakly bound systems (Jones, R.R., You and Bucksbaum 
[1993]). Such pulses can track the wave packet throughout its orbit and detect wave- 
packet motion anywhere in the atom, as was demonstrated experimentally (Raman, 
Conover, Sukenik and Bucksbaum [1996]). In particular, it has been possible to 
monitor the momentum-space probability of the wave packet as a function of time 
(Jones, R.R., [1996]; for details, see also Noordam and Jones, R.R., [1997]). In the 
experimental demonstration of the method also called impulsive momentum retrieval 
Na atoms are first excited by a tunable nanosecond laser to the 25d Rydberg state. 
Subsequently, two half-cycle pulses are applied. Kicking the electron with the first 
pulse, the Rydberg population is redistributed and a complicated, dynamically evolv- 
ing wave packet is created. The second pulse is used to measure the ionization as 
a function of its peak field and delay relative to the first pulse. If the duration of 
the probe pulse is negligible compared with the time scale for variations in the posi- 
tion and momentum of the electron, then the energy gained (or lost) by the electron 
depends only on its initial momentum (i.e., the momentum immediately before the 
application of the pulse) and the time-integrated field. By measuring the ionization 
probability (at threshold) as a function of field, the momentum distribution of the 
initial state along the field axis can be obtained. 

4.5.2. Cyclotron state of a trapped electron 

Let us consider an electron in a Penning trap 125 and assume a uniform magnetic field 
along the positive z axis. The cyclotron and the axial motions are well separated 
in their frequency ranges (for a review, see Brown, L.S., and Gabrielse [1986]). It 
was proposed to reconstruct the cyclotron state by mapping it onto the axial degree 
of freedom (Mancini and Tombesi [1997a]). By applying appropriate fields on the 
trapped electron, a coupling between the two degrees of freedom may be realized 
such that the Hamiltonian in the (cyclotron) interaction picture is of the form 

H = hv z (a\a z + |) + hgx c ((p)z. (238) 

124 For these pulses, commonly referred to as HCP's, see also Greene, Federici, Dykaar, Jones, R.R., 
and Bucksbaum [1991]; You, Jones, R.R., Bucksbaum and Dykaar [1993]. 
125 For the trap, see Penning [1936]. 
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Here, a \ (a z ) is the creation (annihilation) operator for the axial motion of frequency 
u z , and z and x c (ip), respectively, are the axial positon operator and the quadrature 
operator of the cyclotron motion, g being the (field induced) coupling constant. For 
an interaction time that is short compared with the axial period, r <C 2-k/u z , the 
Hamiltonian (238) gives rise to a unitary kick and the time evolution of the mean 
axial momentum yields the mean of the quadrature operator, 



By measuring the current due to the induced charge variation on the cap electrodes 
of the trap, one may obtain the axial momentum and then the cyclotron quadrature 
component. Knowing the cyclotron quadrature-component statistics, the quantum 
state can be derived; e.g., by applying the methods of §§ 3.1., 3.2. and 3.3.1.. 

To determine the cyclotron state from the displaced number statistics (§§ 3.3.2. 
and 3.5.), it was proposed to perform a QND measurement of the cyclotron excita- 
tion number (Mancini and Tombesi [1997a]). Using a magnetic bottle configuration 
(Brown, L.S., and Gabrielse [1986]), an interaction Hamiltonian of the type 



may be realized. It gives rise to a dependence of the axial angular frequency on the 
number of cyclotron excitations, which may be obtained by probing the resonance 
frequency of the output electric signal. The displacement of the cyclotron state may 
be realized by a driving field acting immediately before the measurement process 
induced by the Hamiltonian (240). 

4.5.3. Electron beam 

A method was also proposed for determining the quantum state of the (one-dimens- 
ional) transverse motion of a beam of identically prepared electrons (Tegmark [1996]). 
Originating from the electron source, the beam crosses a magnetic field and in con- 
junction with collimators a highly monochromatic beam is prepared. This beam 
enters a shielded box, wherein a harmonic potential acts in the x-direction. The 
needed potential is obtained by subdividing the walls of the box into a large number 
of metal plates which are insulated from one another, and by fixing the potentials 
of the plates appropriately. A slidable detector (along the z-direction of longitudi- 
nal motion of the beam) allows one to measure the particle density p(x) along the 
x-direction. Performing these measurements for a series of locations of the slidable 
detector, the density matrix characterizing the transverse motion of the electrons of 
the beam may be reconstructed. 

4.6. Spin and angular momentum systems 

In the history of quantum-state measurement spin (angular momentum) systems have 
played a particular role, because of the finite-dimensional Hilbert space (see, e.g., Fano 
[1957] and the references in footnote 41 for the Pauli problem). Spin-quantum-state 
reconstruction is usually based on data obtained from Stern-Gerlach experiments, in 
which particles with magnetic momentum are deflected in a magnetic- field gradient 



{pz(t+r)) = (p z (t)) + hgr(x c (ip)) cos(v z t). 



(239) 



H' = hK,ala c z 2 



(240) 
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and the spin (or angular momentum) projection onto the direction of the field gradient 
is measured (Stern [1921], Gerlach and Stern [1921, 1922]; Feynman, Leighton and 
Sands [1965]; for the optical Stern-Gerlach effect, i.e., state-dependent deflection of 
particles in optical fields, see Kazantsev [1975, 1978], Cook [1978], Tanguy, Reynaud 
and Cohen- Tannoudji [1984], and for its experimental realization, see Sleator, Pfau, 
Balykin, Carnal and Mlynek [1992]). Interacting the particles with fields of magnetic 
multipoles, multipole moments of the particles can be measured (Bohn [1991]). 

In a Stern-Gerlach analyzer as described by Feynman, Leighton and Sands [1965] 
a beam of particles is split in an inhomogeneous magnetic field into different paths 
provided with gates that may be opened or closed. The diagonal density-matrix 
elements g mm are then given by the probability that the particle is deflected into 
the corresponding mth path. Measurement of the off diagonal element Qkn requires 
two analyzers. In the first analyzer, the nth and fcth paths are opened while the 
remaining paths are closed. In the second analyzer the field gradient is in the plane 
that is perpendicular to the quantization axis defined by the first analyzer. Let the 
orientation of the field gradient of the second analyzer be ip. The probability p m ((p) 
of the particle being deflected into the mth path in the second analyzer is then given 
by (Gale, Guth and Trammel [1968]) 

Pm(<f) = 

+2|e fcn |d2(7r/2)t#£(7r/2) cos \{k-n) V + Vkn } }, (241) 

where Qk n = \Qkn\ ex P(^fcn)- The matrix d£2(i?) is the ^-dependent part of the Wigner 
rotation matrix (see Wigner [1959]), and N = (g^ + £> nn ) _1 . Measuring (for known 
Qkk and Q nn ) the probabilities for two different orientations tp, the off diagonal-density 
matrix element Qk n can be determined. 

The quantum state can also be determined by using a single Stern-Gerlach appa- 
ratus with variable orientation (Newton and Young [1968]). The probability p m (i?, (p) 
that the projection of the spin on the quantization axis (■&, ip) is m can be expressed 
in terms of the density matrix elements Qkn as 

Pm(*,<p)= E e-^- fe ^W4 J 2W Qkn, (242) 

n,k=—J 

which can be inverted to obtain 

Qkk+w = E E C m,-m,0 C k+w,-k,w ~JW77T X (^)wm, (243) 
m=-Jj=0 

where 

i 2J 

X($) wm = — — £ ^'PrnfrVs), (244) 
4J + 1 s=-2J 

(p s = 27rs/(4J+ 1), and are Clebsch-Gordan coefficients. It is seen that in 

order to reconstruct the density matrix, it is sufficient to keep t? fixed and take ip at 
4J+1 values such that d^'oW does not vanish for any combination j and w. It should 
be noted that the Stern-Gerlach apparatus can also be allowed to take all possible 



dS(T/2) 



Qkk + 



42(^/2) 
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orientations (Dodonov and Man'ko, V.I., [1997]). In this case one need not search for 
an optimum choice of orientations, since all possible orientations are considered on 
an equal footing. 

As already outlined, the problem of reconstruction of spin density matrices has also 
been treated by means of the principle of maximum entropy (Buzek, Drobny, Adam, 
G., Derka and Knight [1997]; see § 3.9.3.) and Bayesian inference (Jones, K.R.W., 
[1991]; Derka, Buzek and Adam, G., [1996]; Derka, Buzek, Adam, G., and Knight 
[1996]; see § 3.9.4.). It was also proposed to tomographically reconstruct the discrete 
Wigner function of a spin system in analogy to the inverse Radon transform for the 
continuous Wigner function (Leonhardt [1995, 1996]). However, the scheme does not 
only require measurement of probability distributions of spin (angular momentum) 
components but also those of discrete phase states. 

With regard to pure states, it can be shown that the spin wave function can be 
reconstructed from the populations of the spin projections on two directions deviated 
from each other by an infinitesimal angle (Weigert [1992]). Particular attention was 
paid to the determination of the class of pure states which are eigenstates of rotated 
spin projection operators J n = n x J x +n y J y + n z J z , where n= (n x , n y , n z ) is a direction 
vector of unit length. As shown, measurement of spin projections on three direc- 
tions is sufficient to determine the unknown direction n and the corresponding spin 
component eigenvalue (Ivanovic [1993]). 

Finally, experiments for reconstructing the angular momentum density matrix 
of electrons in hydrogen atoms from measurements of the Stark structure of the 
emitted light have been reported (Havener, Rouze, Westerveld and Risley [1986]; 
Ashburn, Cline, Stone, van der Burgt, Westerveld and Risley [1989]; Ashburn, Cline, 
van der Burgt, Westerveld and Risley [1990]; Cline, van der Burgt, Westerveld and 
Risley [1990]; Renwick, Martell, Weaver and Risley [1993]; Seifert, Gibson and Risley 
[1995]). In the experiments, the density matrix is considered for the manifold of 
states of principal quantum number n = 3 and n = 2, the matrix elements being 
parametrized by the orbital and magnetic quantum numbers I and m respectively. 
The excited hydrogen atoms are produced by collisions of protons with noble gas 
atoms. Measuring the coherence parameters of the emitted light [transitions (n = 3) 
— > (^ = 2) and (n = 2) — > {n = 1)] in dependence on the external electric field, one can 
infer the density-matrix elements (for theoretical results, see Jain, Lin and Fritsch 
[1987a,b, 1988]). 

4.7. Crystal lattices 

More than ten years before the first demonstration of reconstruction of the quan- 
tum state of light experiments were already performed for reconstructing off-diagonal 
elements of the single-particle density matrix of crystal lattices (Golovchenko, Ka- 
plan, Kincaid, Levesque, Meixner, Robbins and Felsteiner [1981]; Schiilke, Bonse and 
Mourikis [1981]). The methods are based on Compton scattering or X-ray scattering. 
Before giving the basic idea of these methods, let us outline briefly the information 
available in standard X-ray and Compton scattering experiments. 

In conventional X-ray diffraction the Bragg diffracted intensity corresponding to 
a reciprocal lattice vector g is proportional to the squared modulus of the form factor 
F(g). Provided that the phase of the form factor can be determined, some limited 
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insight into the one-particle density matrix in momentum space, g(p, p + g), can be 
gained. To be more specific, the form factor yields the momentum average of the 
density matrix according to 

F(g) = J d P 0(p,p + g). (245) 

In conventional Compton scattering one may obtain the position average of the den- 
sity matrix in position space, g(r,r'). The one-dimensional Fourier transform of the 
Compton profile, the so-called reciprocal form factor B(r), can be given by 

B(r) = J dr'^r'y + r). (246) 

Obviously, conventional methods do not provide the full density matrix, neither in 
position nor in momentum space. 

The basic idea to overcome this problem consists in the use of coherent scattering 
methods (Golovchenko, Kaplan, Kincaid, Levesque, Meixner, Robbins and Felsteiner 
[1981]; Schiilke, Bonse and Mourikis [1981], Schiilke and Mourikis [1986]; Schiilke 
[1988]). Let us consider the situation in Compton scattering. In conventional experi- 
ments a propagating plane wave is used, so that each point in position space is excited 
equally. This explains the fact that one can only measure position-space averages of 
the density matrix as given in eq. (246). To avoid this limitation, one may introduce 
a weighting of certain positions within the elementary cell. This can be done by using 
a standing-wave field with a spatial periodicity of the nodes that is commensurable 
with the lattice periodicity. The spatial distribution of the corresponding intensity is 
given by 

J(r) = /i + I 2 + 2^hh cos(g-r + </?), (247) 

where the phase tp defines the position of the nodes with respect to the atomic planes of 
the lattice. One may perform two measurements of the Compton profiles by choosing 
the values of the phase to be (p = 0, it. The difference of the corresponding form 
factors, AB(r; g), is then given by 

AB(r; g) ~ J dr' ^(r', r' + r) cos(g-r'). (248) 

That is, coherent Compton scattering allows to determine spatial Fourier transforms 
of the one-particle density matrix. 
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A. Radiation field quantization 

From Maxwell's theory it is well known that the vector potential of the radiation field 
in free space can be given by an expansion in transverse travelling waves. For deal- 
ing with radiation inside resonators, such as cavities bounded by perfectly reflecting 
mirrors, a standing-wave expansion is appropriate, where the spatial structure of the 
waves sensitively depends on the resonator geometry through the boundary condi- 
tions. An expansion of the vector potential 

A(r,t) = A (+) (r,t) + A M (r,t), A«(r,t) = [A (+) (r,t)]* (A.l) 

in orthogonal transverse waves, 

A( + )(r,t) = ^A A (r) e -^a A , (A.2) 

A 

is usually called mode expansion. The mode functions A A (r) satisfy the Helmholtz 
equation 

V x Vx A A -^A A = 0, (A.3) 

cr 

where the condition of transversality, VA A = 0, implies that V x V x A A = — AA A . The 
mode amplitudes a\(t) evolve according to harmonic-oscillator equations of motion, 

a x = -iuj x a x . (A.4) 

Since harmonic oscillators can be associated with the modes, quantization of ra- 
diation reduces to quantization of harmonic oscillators. In eq. (A.l), the complex 
amplitues a A and a* x are regarded as non-Hermitian destruction and creation op- 
erators a A and a\, respectively, which satisfy the well-known (equal-time) bosonic 
commutation relations: 

[a x ,a{,} = S XX ', [a\,ay] = 0. (A. 5) 

The bosonic excitations introduced in this way are called photons. It should be 
pointed out that the expansion in monochromatic modes of the (Heisenberg) operator 
of the vector potential, 

A(r,t) = AW(M) + A H (r,i), (A.6) 
AW(r,t) = ^A A (r)e^a A , (A.7) 

A 

can be used to introduce photons that are associated with other than monochromatic 
waves; viz., 

^ = E^a«a, (A.8) 

A 

A(r, t) = £ [A^r, t) ^ + A>, t) aj] , (A.9) 

where C AtA are elements of a unitary matrix, so that the a'^ and a!^ are again photon 
destruction and creation operators, respectively, and 

A;(r,t) = ]TC M A A (r) e -^ (A.10) 

A 
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(for details, see Titulaer and Glauber [1966]). The spatial-temporal (nonmonochro- 
matic) modes can be used advantageously for describing and analysing pulse-like 
radiation. 

The mirrors used in resonator experiments are of course not perfectly reflecting. 
Moreover, fractionally transparent mirrors are deliberately used to open input and 
output ports in resonator equipments. Another typical example is the use of beam 
splitters in various interference and correlation experiments. In many cases these and 
related (passive) optical instruments can be regarded as macroscopic dielectric bodies 
that respond linearly to radiation and whose action can be included phenomenologi- 
cally in the Maxwell theory through a space-dependent permittivity e (Knoll, Vogel, 
W., and Welsch [1987], Glauber and Lewenstein [1991]). When the radiation under 
consideration is in a frequency intervall where the effect of dispersion is sufficiently 
small, then the dependence of the permittivity on frequency can be disregarded. Per- 
forming a mode expansion (A. 2), the mode functions are now solutions of a modified 
Helmholtz equation (A. 3), with v = v(r) = c/n(r) in place of c, where n(r) = \Je(r) 
is the (real) refractive index. The mode functions satisfy a generalized condition of 
transversality, V-(eA^) =0, and the relation of orthogonality requires inclusion of e 
in the space integral (for details, see Vogel, W., and Welsch [1994]), because of the 
dependence on r of the permittivity. 126 



B. Quantum-state representations 

Among the wide variety of possible radiation-field states, there are some fundamental 
states which, with regard to representation and measurement, play a special role in 
quantum optics. Let us consider a (free) radiation field whose vector potential is 
given by a mode decomposition [eq. (A. 7) or eq. (A. 9)] and restrict our attention to 
single-mode states. The single-mode states can then be used to build up multimode 
states as (direct-)product states. 



B.l. Fock states 

A mode (with photon destruction and creation operators a and a', respectively) is 
said to be in a Fock state (photon-number state) |n), when the state is an eigenstate 
of the photon-number operator n = a^ci, 

n\n)—n\n), n — 0,1,2, (B.l) 

Each number state \n) can be obtained from the vacuum state |0) according to 

H = 4 T ^|0). (B.2) 



i26p or extensions to dispersive and absorbing media, see Gruner and Welsch [1996] and references 
therein. In this case, the permittivity is a complex function of frequency, the real and imaginary 
parts being related to each other according the Kramers-Kronig relations, and a mode decomposition 
outlined here fails. 



B.2. Quadrature-component states 
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Since the photon-number states form an orthogonal and complete set of states, the 
density operator of any (single-mode) quantum state can be expanded in the photon- 
number basis, 

£=$»H£K)(n'|, (B.3) 

n,n' 

the diagonal elements p n — (n\g\n) being the probabilities for observing n photons. 



B.2. Quadrature-component states 

Next let us consider a (single- mode) field-strength operator of the type of F(ip) = 
\F\(ae~ t(p + a)e llf ). Depending on the choice of the (scaling) factor \F\, the operator 
F may represent the vector potential or the electric or magnetic fields associated with 
the mode. In what follows we set F = 2 1 ^ 2 \F\x((f), where 

x(<p) = 2- 1/2 (ae' iLp + tie*) (B.4) 

is called quadrature-component operator. The quadrature components x(ip) and x(ip') 
are observables that cannot be measured simultaneously for (p^(p' + kir, k = 0, ±1, . . ., 
because of [x(ip),x(ip')] 7^0. For chosen ip quadrature-component states \x, ip) can be 
introduced (Schubert and Vogel, W., [1978a,b], see also Vogel, W., and Welsch [1994]), 
solving the eigenvalue problem 

x(ip) \x, ip) — x \x, ip), (B.5) 

—00 < x < 00. Obviously, \x, ip) is related to \x) — \x, if — 0) as \x, ip) = e mip \x), and 
can be given by 



pimp 

\x,<p) = Tr-^e- 2 /^ Hn (x) |n) 

n V2"n! 



7r 1/4 e x2/2 exp 



_i (e i ^) 2 + 2 1 ' 2 xe i ^ 



|0) (B.6) 



[H n (a;), Hermite polynomial]. The states \x,ip) define an orthogonal Hilbert space 
basis for each value of ip, so that the density operator can be expanded as 

g = j dx j dx' \x,ip)(x,ip\g\x',ip)(x',ip\. (B.7) 

In particular, p(x, ip) = (x, tp\g\x, ip) is the probability density for observing the value 
x of the quadrature component x(ip). Note that the symmetry relation p(x, ip + n) = 
p(—x, ip) holds. 



B.3. Coherent states 

The coherent states introduced by Schrodinger [1926] to simulate the motion of a 
(near) classical particle in a harmonic potential are usually defined by the eigenvalue 
problem 127 

a\a) = a \a) (B.8) 

127 The coherent states were recovered in the 1960's (Klauder [1960, 1963a, b]; Glauber [1963a]; 
Sudarshan [1963]; Klauder and Sudarshan [1968]). 
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(a, complex). Equivalently, 

\a) = D(a)\0), (B.9) 
where D(a) is the coherent displacement operator, 

D(a) = exp(aa f - a*a) (B.10) 

[note that D(a)aD\a) — a — a]. The coherent states are not orthogonal, 

(a'\a) = expf — — a'\ 2 ) exp ^(aa'* — a* a') , (B-ll) 



and they are overcomplete (Cahill [1965]; Bacry, Grossmann and Zak [1975]), so that 
there are various possibilities of representing the density operator. In particular, it 
can be expanded as 

g = -^jd 2 aj dV \a)(a\g\a')(a'\ (B.12) 

(d 2 a = d Re a dim a). 

Equation (B.9) can be extented in order to define generalized coherent states by 
the action of the displacement operator on arbitrary states of the Hilbert space (see, 
e.g., Perelomov [1986]). In particular, the action of D(pt) for all values of a on a set 
of states |\1/a) that resolve the unity introduces a foliation of the Hilbert space into 
orbits such that each state belongs exactly to one orbit of generalized coherent states. 
For example, the Fock states can be used to define the generalized coherent states, 

\n,a) = D{a)\n), (B.13) 

which are also called displaced Fock states. 128 Their expansion in the ordinary Fock 
basis reads as 

| re, a) — \m)(m\n, a), (B-14) 

m 

where the expansion coefficients (m\n, a) can be given by 

i { m,n } I I m+n—2l p —i(m—n)(p 

H «, a) = (mtak-W Y £ (-1)- (B.15, 



=0 



[{m, n} = min (m, re), a — \a\e lip ] or 

(m\n,a) = exp (-||«| 2 ) a m ~ n l^- n (\a\ 2 ) (B.16) 



[L™(x), Laguerre polynomial]. 

128 Similar to the coherent states, these states correspond to wave packets of a particle in a harmonic 
potential which keep their shapes and follow the classical motion (Husimi [1953]; Senitzky [1954]; 
Plebanski [1954, 1955, 1956]; Epstein [1959]). 



B.4. s -Parametrized phase-space functions 
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B.4. s-Parametrized phase-space functions 

Representations of the density operator in terms of phase-space functions which are 
formally similar to classical probability distributions are used frequently. Defining an 
s-parameterized displacement operator by 

D(a;s) =e s ^ 2/2 D(a), (B.17) 

the following expansion of the density operator can be proved correct: 



-I' 

7T J 



' ' d 2 aTr{gD{-a;s)}D{a;-s) (B.18) 



(for details, see, e.g., Cahill and Glauber [1969a,b]; Agarwal and Wolf [1970]; Perina 
[1991]). 129 Introducing the Fourier transform of D(a; s), 

5(a -a;s) = \\ d 2 (3D{(3] s) expfa/3* - a*(3\ 

7T J 

= ^) fl(Q '>fcl) D-Ha), (B.19, 
the expansion (B.18) can be rewritten as 130 

g = 7T J d 2 aP(a;s)8(a- a;-s), (B.20) 

where 

P(a; s) = Tr{g5(a -a;s)} (B.21) 

is the s-parametrized phase-space function, which is normalized to unity. 131 Cases of 
particular interest are the Glauber-Sudarshan P function, P(a) = P(a; 1) (Glauber 
[1963b, 1963c]; Sudarshan [1963]), the Wigner function W(a) = P(a;0) (Wigner 
[1932]) and the Q function Q(a) = P{a; -1) (Husimi [1940]). Note that 5(a-a, -1) 
= 7T _1 | a) (a |, which implies that Q(a) = n^ 1 (a\g\a). The expectation value of an 
operator F can then be given by 

(F) = J d 2 a P(a; s) F(a; s), (B.22) 

where F(a; s)/ir is defined according to eq. (B.21) with F in place of g. 
From eqs. (B.19) and (B.21) we see that 

$(ar,s)=Tr{QD(a;s)}, (B.23) 

can be regarded as the characteristic function of the phase-space function P(a; s), 

P(a; s) = J d 2 p exp(a/T - a*/3)$(/3; s). (B.24) 

i29p or s = o the expansion (B.18) is known from Weyl's quantization method (Weyl [1927]). 

130 Note that eqs. (B.20) and (B.21) also apply to other than density operators. Let F be an 
operator function of a and a'. Then it can be shown that F — J d 2 a F(a;s) 5(a — a;s), where 
F(a; s) =7rTr{F 8 (a — a; — s)} is the associated c-number function in s order. 

131 Frequently the definition P{q,p: s) = 2~ 1 P[a = 2~ 1 / 2 (q + ip)] is used, so that P(q,p; s) can be 
regarded as a function of "position" and "momentum" , with J dq J dpP(q,p; s) = 1. 
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Using eqs. (B.17), (B.23), and (B.24), the phase-space functions P(a; s) can be related 
to each other as 

P(a; s) = ^J d 2 {3 exp[l(s-s')\(3\ 2 

x J d 2 -fexp[(a--f)(3*-(a*- 1 *)f3]P(r,s'). (B.25) 

It should be noted that the formalism can be extended to finite dimensional Hilbert 
spaces (Opatrny, Welsch and Buzek [1996]; for the Wigner function, see Wootters 
[1987], Leonhardt [1995,1996]; for the Q function, see Opatrny, Buzek, Bajer and 
Drobny [1995] and Galetti and Marchiolli [1996]). 132 

B.5. Quantum state and quadrature components 

It is worth noting that there is a one-to-one correspondence between the phase-space 
function P(a; s) and the set of quadrature-component distributions p(x, tp) for all 
values of tp within a tt interval (Vogel, K., and Risken [1989]). Introducing the char- 
acteristic function ty(z,tp) of the quadrature-component distribution p(x,tp), 

V(z,tp) = J dxe izx p(x,tp), (B.26) 

it can be shown that 

*(z,<p) = e^ 2/4 $(^2- 1/2 e^;s). (B.27) 
Equations (B.24), (B.26), and (B.27) reveal that 

1 t 

p(x,tp) — — J dz exp(^—izx — \sz 2 ^j 

x J d 2 a exp[iz(a\x((p)\a)]P(a;s) (B.28) 

and 

x exp[-iz{a\x(ip)\a)] J dxe tzx p(x, </?)|. (B.29) 

In other words, knowledge of p(x, ip) for all values of tp within a n interval is equivalent 
to knowledge of the quantum state. 133 In particular, the expansion of the density 
operator as given in eq. (B.18) can be rewritten, on using eqs. (B.23) and (B.27), as 

g = — j % dtp I dz\z\e sz2/ ^(-z,p)]D{iz2- 1 / 2 e^;-s) , (B.30) 
2,71 Jo J V ' 

which offers the possibility of relating g to p(x, tp) as 

q — j dp J dx K(x, tp) p(x, tp), (B.31) 

132 For an approach to phase-space functions for systems with finite-dimensional Hilbert spaces 
which uses a continuous phase space, see Varilly and Gracia-Bondfa [1989]; Dowling, Agarwal and 
Schleich [1994]. 

133 Note that the tp integral in cq. (B.29) can be performed over any n interval. 
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provided that the operator kernel 

K(x,ip) — — ( dz\z\ exp{iz [x(<p) —x]}, (B.32) 
exists (D'Ariano [1995]; D'Ariano, Leonhardt and Paul, H., [1995]). 

C. Photodetection 

In order to give an introduction into the quantum theory of photoelectric detection of 
light (Mandel [1958, 1963], Kelley and Kleiner [1964], Glauber [1965]), let us consider, 
as a simple example of a photodetection device, a large sample of atomic systems that 
are capable of absorbing light through the photoemission of electrons in a certain time 
interval t, t + At. Next suppose that the photoemission of the actual number m of 
photoelectrons is dominated by the process in which exactly m atomic systems are 
involved, so that each emits just one electron. We further assume that the total num- 
ber N of atomic systems is much larger compared with the mean number of emitted 
electrons, so that for any (relevant) actual value m the inequality m<N may be as- 
sumed. Under these assumptions, the main features of the theory may be developed 
by applying Dirac's perturbation theory to the basic process of light absorption and 
combining the corresponding results with methods of classical statistics with respect 
to the ensemble of photoelectrons generated by the absorption processes (see, e.g., 
Vogel, W., and Welsch [1994]). 

Here, we will give a more intuitive analysis rather than an exact derivation. When 
each photon that falls - in the chosen time interval - on the detector gives rise to 
exactly one emitted electron, then the number of counted electrons agrees exactly 
with the number of photons, and the statistics of the counted electrons reflect exactly 
the photon- number statistics; i.e., the probability P rn of detecting m photoelectrons 
is equal to the probability p m of m photons being in the field, 

P m = Pm = (m\g\m). (C.l) 

However owing to losses, the probability r) of converting a photon to an electron is 
less than unity in general (0 < rj < 1). This probability is also called detection (or 
quantum) efficiency Since (under the assumptions made) the individual events of 
emission of a photoelectron can be regarded as being independent of each other, the 
probability P m \ n {rj) of observing m photoelectrons under the condition that n photons 
are present corresponds to a Bernoulli process, 

Pm\n(v)= ( n )v m (l-v) n ' m if rn<n, (C.2) 

and P m \ n {v) = if m > n. The joint probability that n photons are present and 
m photoelectrons are counted is then P m \n(v)Pn, and hence the prior probability of 
detecting m photoelectrons is given by 



00 'n 



P m = J2 P rn\n(v)Pn= £ [ ~ VT^Pn (C.3) 



n=m 
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which for rj — > 1 (perfect detection) reduces to eq. (C.l). Note that eq. (C.3) can be 
inverted in order to obtain p n from P r 

-1,0 

Equation (C.3) can be rewritten as 134 
and the characteristic function 

n(y) 

can be given by 

:), (C.7) 

where the symbol : : indicates normal ordering, with the operator a) to the left of 
the operator a. So far, we have assumed that the photons are effectively associated 
with a single (nonmonochromatic) mode. The extension of the above given formulae 
to multimode fields is straightforward. In particular, the multi-mode version of the 
single-mode photocounting formula (C.5) reads as 

Prn^-lU--^^^), (C.8) 
\ k " Lk - I 

where n fc is the photon-number operator of the fcth mode, and r/ k is the quantum 
efficiency which with the photons of this mode are detected. 



by simply replacing rj with rj 




(C.4) 



(C.5) 



= Y / P m e mv (C.6) 



: exp 



rjh(e iy - 1) 



D. Elements of least-squares inversion 35 

Let f be a (possibly unknown) no dimensional "state" vector and consider a stochastic 
linear transform, 

y = Af + n, (D.l) 

yielding an mo dimensional (mo > no) "data" vector y available from measurements. 
Here, A is a given mo x no matrix, and n is an mo dimensional vector whose random 
elements with zero means and covariance matrix W~ l describe the noise associated 
with realistic measurements. The probability of a state vector f being realized under 
the condition that there is a data vector y can be given by 

P(f|y)~P(y|f)P(f), (D.2) 

134 Note that the idendity \n)(n\ — : (n n /n\)e~ n : is valid, which can easily be proved correct from 
the associated c— number function |(a|n)| 2 of \n)(n\ in normal order. 

135 For details, see, e.g., Golub and van Loan [1989]; Robinson [1991]; Press, Teukolsky, Vetterling 
and Flannery [1995]. 
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where for Gaussian noise the conditional probability P(y\f) of observing the data 
vector y is 

P(y|f) ~ exp[-i(y - A^W(y - Af)] . (D.3) 

The probability P(f) is a measure of a priori knowledge of the state and can be set 
constant if no state is preferred. The most probable state f can then be found by 
minimization of 136 

C(f) = (y-Af)t^(y-Af), (D.4) 

from which 

A^WAi = A^Wy, (D.5) 
and when A^VFA is not singular, then 

f = {A^W Ay 1 A^Wy . (D.6) 

Otherwise, the inversion of eq. (D.5) is not unique and further criteria must be used 
to select a solution. When the matrix W is not known, then it may be set a multiple 
of a unity matrix, so that eq. (D.6) reduces to 

f = (AUr^V, (D.7) 

provided that A^A is not singular. Equation (D.7) still gives the correct averaged 
inversion, but the statistical fluctuation of the result may be (slightly) enhanced. 

If the data are not sensitive enough to some state-vector components, then these 
components can hardly be determined with reasonable accuracy. Mathematically, 
i^WA becomes (quasi-) singular and regularizations, such as Tikhonov regulariza- 
tion and singular- value decomposition, are required to solve approximately eq. (D.5). 
For simplicity let us set W = I (I, unity matrix). Using Tikhonov regularization, it 
is assumed that some components of the state vector can be preferred by a properly 
chosen a priori probability P(f), such as 

P(f) ~ exp(-iA 2 ftf) , (D.8) 

the parameter A (A>0) being a measure of the strength of regularization. Maximiza- 
tion of P(f |y) then yields, on recalling eqs. (D.2) and (D.3), 

f =(A 2 I + i4 t A)" 1 A t y- (D.9) 

Note that A 2 I + A* A has only positive eigenvalues and is thus always invertible. A 
possible choice of A is based on the so-called L curve, which is a log-log plot of ||f|| 
versus ||Ay||, Ay = y — Af, for different values of A. The points on the horizontal 
branch correspond to large noise, whereas the points on the vertical branch correspond 
to large data misfit. Optimum choice of A corresponds to points near the corner of 
the L curve. 



!36jf - 1S diagonal (i.e., the noise is uncorrelated), then eq. (D.4) represents a sum of weighted 
squares of the differences between the components of the data vector y and the components of the 
transformed vector Af, each term of the sum being multiplied by a weight given by the corresponding 
element of W. 
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Applying singular- value decomposition, the inversion of the matrix A* A is per- 
formed such that their eigenvalues whose absolute values are smaller than the (pos- 
itive) parameter of regularization <t are treated as zeros, but the inversions are set 
zero (instead to infinity). This operation is called "pseudoinverse" of a matrix, 

f = Pseudoinverse(A f A; a )A ] y. (D.10) 

For <7 close to zero the result of eq. (D.10) is similar to that of eq. (D.7). With 
increasing o"o, smaller absolute values of components of f are preferred. 

The effect of the regularization parameters A and a is similar. The statistical error 
of the reconstructed state vector f is decreased, but bias towards zero is produced 
simultaneously. Hence, optimum parameters are those for which the bias is just below 
the statistical fluctuation. The bias can be estimated, e.g., by Monte Carlo generating 
new sets of "synthetic" data from the reconstructed state. From these sets one can 
again reconstruct new sets of f . The difference between the mean value of the states 
reconstructed from the synthetic data and the originally reconstructed state estimates 
the bias. 
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